package org.lestr.secuel.sogl.maths;

import org.lestr.secuel.m3ds.maths.IVecteur3d;

public class Transformation3D {

    double[] mat = new double[16];
    //double[] rot = new double[9];
    //double[] scales = new double[3];
    // allocate the memory only when it is needed. Following three places will allocate the memory,
    // void setScaleTranslation(), void computeScales() and void computeScaleRotation()
    double[] rot = null;
    double[] scales = null;

    // Unknown until lazy classification is done
    private int type = 0;

    // Dirty bit for classification, this is used
    // for classify()
    private static final int AFFINE_BIT     = 0x01;
    private static final int ORTHO_BIT = 0x02;
    private static final int CONGRUENT_BIT  = 0x04;
    private static final int RIGID_BIT      = 0x08;
    private static final int CLASSIFY_BIT   = 0x10;

    // this is used for scales[], rot[]
    private static final int SCALE_BIT      = 0x20;
    private static final int ROTATION_BIT   = 0x40;
    // set when SVD renormalization is necessary
    private static final int SVD_BIT        = 0x80;

    private static final int CLASSIFY_ALL_DIRTY = AFFINE_BIT |
                                                  ORTHO_BIT |
                                                  CONGRUENT_BIT |
                                                  RIGID_BIT |
                                                  CLASSIFY_BIT;
    private static final int ROTSCALESVD_DIRTY = SCALE_BIT |
                                                  ROTATION_BIT |
                                                  SVD_BIT;
    private static final int ALL_DIRTY = CLASSIFY_ALL_DIRTY | ROTSCALESVD_DIRTY;

    private int dirtyBits;

    boolean autoNormalize = false;	// Don't auto normalize by default
    /*
    // reused temporaries for compute_svd
    private boolean svdAllocd =false;
    private double[] u1 = null;
    private double[] v1 = null;
    private double[] t1 = null; // used by both compute_svd and compute_qr
    private double[] t2 = null; // used by both compute_svd and compute_qr
    private double[] ts = null;
    private double[] svdTmp = null;
    private double[] svdRot = null;
    private double[] single_values = null;
    private double[] e = null;
    private double[] svdScales = null;
    // from svrReorder
    private int[] svdOut = null;
    private double[] svdMag = null;

    // from compute_qr
    private double[]   cosl  = null;
    private double[]   cosr  = null;
    private double[]   sinl  = null;
    private double[]   sinr  = null;
    private double[]   qr_m  = null;
    */

    private static final double EPS = 1.110223024E-16;

    static final double EPSILON = 1.0e-10;
    static final double EPSILON_ABSOLUTE = 1.0e-5;
    static final double EPSILON_RELATIVE = 1.0e-4;
    /**
     * A zero matrix.
     */
    public static final int ZERO = 0x01;

   /**
    * An identity matrix.
    */
    public static final int IDENTITY = 0x02;


   /**
    * A Uniform scale matrix with no translation or other
    * off-diagonal components.
    */
    public static final int SCALE = 0x04;

   /**
    * A translation-only matrix with ones on the diagonal.
    *
    */
    public static final int TRANSLATION = 0x08;

   /**
    * The four row vectors that make up an orthogonal matrix form a basis,
    * meaning that they are mutually orthogonal; an orthogonal matrix with
    * positive determinant is a pure rotation matrix; a negative
    * determinant indicates a rotation and a reflection.
    */
    public static final int ORTHOGONAL = 0x10;

   /**
    * This matrix is a rotation and a translation with unity scale;
    * The upper 3x3 of the matrix is orthogonal, and there is a
    * translation component.
    */
    public static final int RIGID = 0x20;

   /**
    * This is an angle and length preserving matrix, meaning that it
    * can translate, rotate, and reflect
    * about an axis, and scale by an amount that is uniform in all directions.
    * These operations preserve the distance between any two points and the
    * angle between any two intersecting lines.
    */
    public static final int CONGRUENT = 0x40;

   /**
    * An affine matrix can translate, rotate, reflect, scale anisotropically,
    * and shear.  Lines remain straight, and parallel lines remain parallel,
    * but the angle between intersecting lines can change. In order for a
    * transform to be classified as affine, the 4th row must be: [0, 0, 0, 1].
    */
    public static final int AFFINE = 0x80;

   /**
    * This matrix has a negative determinant; an orthogonal matrix with
    * a positive determinant is a rotation matrix; an orthogonal matrix
    * with a negative determinant is a reflection and rotation matrix.
    */
    public static final int NEGATIVE_DETERMINANT = 0x100;

    /**
     * The upper 3x3 column vectors that make up an orthogonal
     * matrix form a basis meaning that they are mutually orthogonal.
     * It can have non-uniform or zero x/y/z scale as long as
     * the dot product of any two column is zero.
     * This one is used by Java3D internal only and should not
     * expose to the user.
     */
    private static final int ORTHO = 0x40000000;

    /**
     * Constructs and initializes a transform to the identity matrix.
     */
    public Transformation3D() {
        setIdentity();			// this will also classify the matrix
    }

   /**
     * Returns the type of this matrix as an or'ed bitmask of
     * of all of the type classifications to which it belongs.
     * @return  or'ed bitmask of all of the type classifications
     * of this transform
     */
    public final int getType() {
	if ((dirtyBits & CLASSIFY_BIT) != 0) {
	    classify();
	}
	// clear ORTHO bit which only use internally
	return (type & ~ORTHO);
    }

    // True if type is ORTHO
    // Since ORTHO didn't take into account the last row.
    final boolean isOrtho() {
	if ((dirtyBits & ORTHO_BIT) != 0) {
	    if ((almostZero(mat[0]*mat[2] + mat[4]*mat[6] +
			    mat[8]*mat[10]) &&
		 almostZero(mat[0]*mat[1] + mat[4]*mat[5] +
			    mat[8]*mat[9]) &&
		 almostZero(mat[1]*mat[2] + mat[5]*mat[6] +
			    mat[9]*mat[10]))) {
		type |= ORTHO;
		dirtyBits &= ~ORTHO_BIT;
		return true;
	    } else {
		type &= ~ORTHO;
		dirtyBits &= ~ORTHO_BIT;
		return false;
	    }
	}
	return ((type & ORTHO) != 0);
    }

    final boolean isCongruent() {
	if ((dirtyBits & CONGRUENT_BIT) != 0) {
	    // This will also classify AFFINE
		classifyRigid();
	}
	return ((type & CONGRUENT) != 0);
    }

    final boolean isAffine() {
	if ((dirtyBits & AFFINE_BIT) != 0) {
	    classifyAffine();
	}
	return ((type & AFFINE) != 0);
    }

    final boolean isRigid() {
	if ((dirtyBits & RIGID_BIT) != 0) {


	    // This will also classify AFFINE & CONGRUENT
	    if ((dirtyBits & CONGRUENT_BIT) != 0) {
		classifyRigid();
	    } else {

		if ((type & CONGRUENT) != 0) {
		    // Matrix is Congruent, need only
		    // to check scale
		    double s;
		    if ((dirtyBits & SCALE_BIT) != 0){
			s = mat[0]*mat[0] + mat[4]*mat[4] +
			    mat[8]*mat[8];
			// Note that
			// scales[0] = sqrt(s);
			// but since sqrt(1) = 1,
			// we don't need to do s = sqrt(s) here.
		    } else {
			if(scales == null)
			    scales = new double[3];
			s = scales[0];
		    }
		    if (almostOne(s)) {
			type |= RIGID;
		    } else {
			type &= ~RIGID;
		    }
		} else {
		    // Not even congruent, so isRigid must be false
		    type &= ~RIGID;
		}
		dirtyBits &= ~RIGID_BIT;
	    }
	}
	return ((type & RIGID) != 0);
    }


   /**
     * Returns the least general type of this matrix; the order of
     * generality from least to most is: ZERO, IDENTITY,
     * SCALE/TRANSLATION, ORTHOGONAL, RIGID, CONGRUENT, AFFINE.
     * If the matrix is ORTHOGONAL, calling the method
     * getDeterminantSign() will yield more information.
     * @return the least general matrix type
     */
    public final int getBestType() {
	getType();   // force classify if necessary

	if ((type & ZERO)                 != 0 ) return ZERO;
	if ((type & IDENTITY)             != 0 ) return IDENTITY;
	if ((type & SCALE)                != 0 ) return SCALE;
	if ((type & TRANSLATION)          != 0 ) return TRANSLATION;
	if ((type & ORTHOGONAL)           != 0 ) return ORTHOGONAL;
	if ((type & RIGID)                != 0 ) return RIGID;
	if ((type & CONGRUENT)            != 0 ) return CONGRUENT;
	if ((type & AFFINE)               != 0 ) return AFFINE;
	if ((type & NEGATIVE_DETERMINANT) != 0 ) return NEGATIVE_DETERMINANT;
	return 0;
    }

    /*
    private void print_type() {
        if ((type & ZERO)                 > 0 ) System.err.print(" ZERO");
	if ((type & IDENTITY)             > 0 ) System.err.print(" IDENTITY");
	if ((type & SCALE)                > 0 ) System.err.print(" SCALE");
	if ((type & TRANSLATION)          > 0 ) System.err.print(" TRANSLATION");
	if ((type & ORTHOGONAL)           > 0 ) System.err.print(" ORTHOGONAL");
	if ((type & RIGID)                > 0 ) System.err.print(" RIGID");
	if ((type & CONGRUENT)            > 0 ) System.err.print(" CONGRUENT");
	if ((type & AFFINE)               > 0 ) System.err.print(" AFFINE");
	if ((type & NEGATIVE_DETERMINANT) > 0 ) System.err.print(" NEGATIVE_DETERMINANT");
	}
    */

    /**
     * Returns the sign of the determinant of this matrix; a return value
     * of true indicates a non-negative determinant; a return value of false
     * indicates a negative determinant. A value of true will be returned if
     * the determinant is NaN. In general, an orthogonal matrix
     * with a positive determinant is a pure rotation matrix; an orthogonal
     * matrix with a negative determinant is a both a rotation and a
     * reflection matrix.
     * @return  determinant sign : true means non-negative, false means negative
     */
    public final boolean getDeterminantSign() {
        double det = determinant();
        if (Double.isNaN(det)) {
            return true;
        }
        return det >= 0;
    }

    /**
     * Sets a flag that enables or disables automatic SVD
     * normalization.  If this flag is enabled, an automatic SVD
     * normalization of the rotational components (upper 3x3) of this
     * matrix is done after every subsequent matrix operation that
     * modifies this matrix.  This is functionally equivalent to
     * calling normalize() after every subsequent call, but may be
     * less computationally expensive.
     * The default value for this parameter is false.
     * @param autoNormalize  the boolean state of auto normalization
     */
    public final void setAutoNormalize(boolean autoNormalize) {
	this.autoNormalize = autoNormalize;

	if (autoNormalize) {
	    normalize();
	}
    }

    /**
     * Returns the state of auto-normalization.
     * @return  boolean state of auto-normalization
     */
    public final boolean getAutoNormalize() {
	return this.autoNormalize;
    }
   
    
    
    private static final boolean almostZero(double a) {
	return ((a < EPSILON_ABSOLUTE) && (a > -EPSILON_ABSOLUTE));
    }

    private static final boolean almostOne(double a) {
	return ((a < 1+EPSILON_ABSOLUTE) && (a > 1-EPSILON_ABSOLUTE));
    }

    private static final boolean almostEqual(double a, double b) {
	double diff = a-b;

	if (diff >= 0) {
	    if (diff < EPSILON) {
		return true;
	    }
	    // a > b
	    if ((b > 0) || (a > -b)) {
		return (diff < EPSILON_RELATIVE*a);
	    } else {
		return (diff < -EPSILON_RELATIVE*b);
	    }

	} else {
	    if (diff > -EPSILON) {
		return true;
	    }
	    // a < b
	    if ((b < 0) || (-a > b)) {
		return (diff > EPSILON_RELATIVE*a);
	    } else {
		return (diff > -EPSILON_RELATIVE*b);
	    }
	}
    }

    private final void classifyAffine() {
        if (!isInfOrNaN() &&
                almostZero(mat[12]) &&
                almostZero(mat[13]) &&
                almostZero(mat[14]) &&
                almostOne(mat[15])) {
	    type |= AFFINE;
	} else {
	    type &= ~AFFINE;
	}
	dirtyBits &= ~AFFINE_BIT;
    }

    // same amount of work to classify rigid and congruent
    private final void classifyRigid() {

	if ((dirtyBits & AFFINE_BIT) != 0) {
	    // should not touch ORTHO bit
	    type &= ORTHO;
	    classifyAffine();
	} else {
	    // keep the affine bit if there is one
	    // and clear the others (CONGRUENT/RIGID) bit
	    type &= (ORTHO | AFFINE);
	}

	if ((type & AFFINE) != 0) {
	    // checking orthogonal condition
	    if (isOrtho()) {
		if ((dirtyBits & SCALE_BIT) != 0) {
		    double s0 = mat[0]*mat[0] + mat[4]*mat[4] +
			mat[8]*mat[8];
		    double s1 = mat[1]*mat[1] + mat[5]*mat[5] +
			mat[9]*mat[9];
		    if (almostEqual(s0, s1)) {
			double s2 = mat[2]*mat[2] + mat[6]*mat[6] +
			    mat[10]*mat[10];
			if (almostEqual(s2, s0)) {
			    type |= CONGRUENT;
			    // Note that scales[0] = sqrt(s0);
			    if (almostOne(s0)) {
				type |= RIGID;
			    }
			}
		    }
		} else {
		    if(scales == null)
			scales = new double[3];

		    double s = scales[0];
		    if (almostEqual(s, scales[1]) &&
			almostEqual(s, scales[2])) {
			type |= CONGRUENT;
			if (almostOne(s)) {
			    type |= RIGID;
			}
		    }
		}
	    }
	}
	dirtyBits &= (~RIGID_BIT | ~CONGRUENT_BIT);
    }

    /**
     * Classifies a matrix.
     */
    private final void classify() {

	if ((dirtyBits & (RIGID_BIT|AFFINE_BIT|CONGRUENT_BIT)) != 0) {
	    // Test for RIGID, CONGRUENT, AFFINE.
	    classifyRigid();
	}

	// Test for ZERO, IDENTITY, SCALE, TRANSLATION,
	// ORTHOGONAL, NEGATIVE_DETERMINANT
	if ((type & AFFINE) != 0) {
	    if ((type & CONGRUENT) != 0) {
		if ((type & RIGID) != 0) {
		    if (zeroTranslation()) {
			type |= ORTHOGONAL;
			if (rotateZero()) {
			    // mat[0], mat[5], mat[10] can be only be
			    // 1 or -1 when reach here
			    if ((mat[0] > 0) &&
				(mat[5] > 0) &&
				(mat[10] > 0)) {
				type |= IDENTITY|SCALE|TRANSLATION;
			    }
			}
		    } else {
			if (rotateZero()) {
			    type |= TRANSLATION;
			}
		    }
		} else {
		    // uniform scale
		    if (zeroTranslation() && rotateZero()) {
			type |= SCALE;
		    }
		}

	    }
	} else {
	    // last row is not (0, 0, 0, 1)
	    if (almostZero(mat[12]) &&
		almostZero(mat[13]) &&
		almostZero(mat[14]) &&
		almostZero(mat[15]) &&
		zeroTranslation() &&
		rotateZero() &&
		almostZero(mat[0]) &&
		almostZero(mat[5]) &&
		almostZero(mat[10])) {
		type |= ZERO;
	    }
	}

	if (!getDeterminantSign()) {
	    type |= NEGATIVE_DETERMINANT;
	}
	dirtyBits &= ~CLASSIFY_BIT;
    }

    final boolean zeroTranslation() {
	return (almostZero(mat[3]) &&
		almostZero(mat[7]) &&
		almostZero(mat[11]));
    }

    final boolean rotateZero() {
	return (almostZero(mat[1]) && almostZero(mat[2]) &&
		almostZero(mat[4]) && almostZero(mat[6]) &&
		almostZero(mat[8]) && almostZero(mat[9]));
    }

   /**
     * Returns the matrix elements of this transform as a string.
     * @return  the matrix elements of this transform
     */
    @Override
	public String toString() {
	// also, print classification?
	return
	    mat[0] + ", " + mat[1] + ", " + mat[2] + ", " + mat[3] + "\n" +
	    mat[4] + ", " + mat[5] + ", " + mat[6] + ", " + mat[7] + "\n" +
	    mat[8] + ", " + mat[9] + ", " + mat[10] + ", " + mat[11] + "\n" +
	    mat[12] + ", " + mat[13] + ", " + mat[14] + ", " + mat[15]
	    + "\n";
    }

    /**
     * Sets this transform to the identity matrix.
     */
    public final void setIdentity() {
	mat[0] = 1.0;  mat[1] = 0.0;  mat[2] = 0.0;  mat[3] = 0.0;
	mat[4] = 0.0;  mat[5] = 1.0;  mat[6] = 0.0;  mat[7] = 0.0;
	mat[8] = 0.0;  mat[9] = 0.0;  mat[10] = 1.0; mat[11] = 0.0;
	mat[12] = 0.0; mat[13] = 0.0; mat[14] = 0.0; mat[15] = 1.0;
	type = IDENTITY | SCALE |  ORTHOGONAL | RIGID | CONGRUENT |
	       AFFINE | TRANSLATION | ORTHO;
	dirtyBits = SCALE_BIT | ROTATION_BIT;
	// No need to set SVD_BIT
    }

   /**
     * Sets this transform to all zeros.
     */
    public final void setZero() {
	mat[0] = 0.0;  mat[1] = 0.0;  mat[2] = 0.0;  mat[3] = 0.0;
	mat[4] = 0.0;  mat[5] = 0.0;  mat[6] = 0.0;  mat[7] = 0.0;
	mat[8] = 0.0;  mat[9] = 0.0;  mat[10] = 0.0; mat[11] = 0.0;
	mat[12] = 0.0; mat[13] = 0.0; mat[14] = 0.0; mat[15] = 0.0;

	type = ZERO | ORTHO;
	dirtyBits = SCALE_BIT | ROTATION_BIT;
    }


   /**
     * Adds this transform to transform t1 and places the result into
     * this: this = this + t1.
     * @param t1  the transform to be added to this transform
     */
    public final void add(Transformation3D t1) {
	for (int i=0 ; i<16 ; i++) {
	    mat[i] += t1.mat[i];
	}

	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}

    }

    /**
     * Adds transforms t1 and t2 and places the result into this transform.
     * @param t1  the transform to be added
     * @param t2  the transform to be added
     */
    public final void add(Transformation3D t1, Transformation3D t2) {
	for (int i=0 ; i<16 ; i++) {
	    mat[i] = t1.mat[i] + t2.mat[i];
	}

	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}

    }

    /**
     * Subtracts transform t1 from this transform and places the result
     * into this: this = this - t1.
     * @param t1  the transform to be subtracted from this transform
     */
    public final void sub(Transformation3D t1) {
	for (int i=0 ; i<16 ; i++) {
	    mat[i] -= t1.mat[i];
	}

	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}
    }


    /**
     * Subtracts transform t2 from transform t1 and places the result into
     * this: this = t1 - t2.
     * @param t1   the left transform
     * @param t2   the right transform
     */
    public final void sub(Transformation3D t1, Transformation3D t2) {
	for (int i=0 ; i<16 ; i++) {
	    mat[i] = t1.mat[i] - t2.mat[i];
	}

	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}

    }


   /**
     * Transposes this matrix in place.
     */
    public final void transpose() {
        double temp;

        temp = mat[4];
        mat[4] = mat[1];
        mat[1] = temp;

        temp = mat[8];
        mat[8] = mat[2];
        mat[2] = temp;

        temp = mat[12];
        mat[12] = mat[3];
        mat[3] = temp;

        temp = mat[9];
        mat[9] = mat[6];
        mat[6] = temp;

        temp = mat[13];
        mat[13] = mat[7];
        mat[7] = temp;

        temp = mat[14];
        mat[14] = mat[11];
        mat[11] = temp;

	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}
    }

    /**
     * Transposes transform t1 and places the value into this transform.
     * The transform t1 is not modified.
     * @param t1  the transform whose transpose is placed into this transform
     */
    public final void transpose(Transformation3D t1) {

       if (this != t1) {
           mat[0] =  t1.mat[0];
           mat[1] =  t1.mat[4];
           mat[2] =  t1.mat[8];
           mat[3] =  t1.mat[12];
           mat[4] =  t1.mat[1];
           mat[5] =  t1.mat[5];
           mat[6] =  t1.mat[9];
           mat[7] =  t1.mat[13];
           mat[8] =  t1.mat[2];
           mat[9] =  t1.mat[6];
           mat[10] = t1.mat[10];
           mat[11] = t1.mat[14];
           mat[12] = t1.mat[3];
           mat[13] = t1.mat[7];
           mat[14] = t1.mat[11];
           mat[15] = t1.mat[15];

	   dirtyBits = ALL_DIRTY;

	   if (autoNormalize) {
	       normalize();
	   }
       } else {
           this.transpose();
       }

    }

   /**
     * Sets the rotational component (upper 3x3) of this transform to the
     * matrix equivalent values of the axis-angle argument; the other
     * elements of this transform are unchanged; any pre-existing scale
     * in the transform is preserved.
     * @param a1 the axis-angle to be converted (x, y, z, angle)
     */
    public final void setRotation(IVecteur3d u, double angle) {

	if ((dirtyBits & SCALE_BIT)!= 0) {
	    computeScales(false);
	}

	double mag = Math.sqrt(u.obtenirX()*u.obtenirX() + u.obtenirY()*u.obtenirY() + u.obtenirZ()*u.obtenirZ());

	if (almostZero(mag)) {
	    mat[0] = scales[0];
	    mat[1] = 0.0;
	    mat[2] = 0.0;
	    mat[4] = 0.0;
	    mat[5] = scales[1];
	    mat[6] = 0.0;
	    mat[8] = 0.0;
	    mat[9] = 0.0;
	    mat[10] = scales[2];
	} else {
	    mag = 1.0/mag;
	    double ax = u.obtenirX()*mag;
	    double ay = u.obtenirY()*mag;
	    double az = u.obtenirZ()*mag;

	    double sinTheta = Math.sin(angle);
	    double cosTheta = Math.cos(angle);
	    double t = 1.0 - cosTheta;

	    double xz = ax * az;
	    double xy = ax * ay;
	    double yz = ay * az;

	    mat[0] = (t * ax * ax + cosTheta)*scales[0];
	    mat[1] = (t * xy - sinTheta * az)*scales[1];
	    mat[2] = (t * xz + sinTheta * ay)*scales[2];

	    mat[4] = (t * xy + sinTheta * az)*scales[0];
	    mat[5] = (t * ay * ay + cosTheta)*scales[1];
	    mat[6] = (t * yz - sinTheta * ax)*scales[2];

	    mat[8] = (t * xz - sinTheta * ay)*scales[0];
	    mat[9] = (t * yz + sinTheta * ax)*scales[1];
	    mat[10] = (t * az * az + cosTheta)*scales[2];
	}

        // Issue 253: set all dirty bits if input is infinity or NaN
        if (isInfOrNaN(u)) {
            dirtyBits = ALL_DIRTY;
            return;
        }

	// Rigid remain rigid, congruent remain congruent after
	// set rotation
	dirtyBits |= CLASSIFY_BIT | ROTATION_BIT;
	dirtyBits &= ~ORTHO_BIT;
	type |= ORTHO;
	type &= ~(ORTHOGONAL|IDENTITY|SCALE|TRANSLATION|SCALE|ZERO);
    }


    /**
     * Sets the value of this transform to a counter clockwise rotation
     * about the x axis. All of the non-rotational components are set as
     * if this were an identity matrix.
     * @param angle the angle to rotate about the X axis in radians
     */
    public void rotX(double angle) {
        double sinAngle = Math.sin(angle);
        double cosAngle = Math.cos(angle);

        mat[0] = 1.0;
        mat[1] = 0.0;
        mat[2] = 0.0;
        mat[3] = 0.0;

        mat[4] = 0.0;
        mat[5] = cosAngle;
        mat[6] = -sinAngle;
        mat[7] = 0.0;

        mat[8] = 0.0;
        mat[9] = sinAngle;
        mat[10] = cosAngle;
        mat[11] = 0.0;

        mat[12] = 0.0;
        mat[13] = 0.0;
        mat[14] = 0.0;
        mat[15] = 1.0;

        // Issue 253: set all dirty bits if input is infinity or NaN
        if (isInfOrNaN(angle)) {
            dirtyBits = ALL_DIRTY;
            return;
        }

	type = CONGRUENT | AFFINE | RIGID | ORTHO;
	dirtyBits = CLASSIFY_BIT | ROTATION_BIT | SCALE_BIT;
    }

    /**
     * Sets the value of this transform to a counter clockwise rotation about
     * the y axis. All of the non-rotational components are set as if this
     * were an identity matrix.
     * @param angle the angle to rotate about the Y axis in radians
     */
    public void rotY(double angle) {
        double sinAngle = Math.sin(angle);
        double cosAngle = Math.cos(angle);

        mat[0] = cosAngle;
        mat[1] = 0.0;
        mat[2] = sinAngle;
        mat[3] = 0.0;

        mat[4] = 0.0;
        mat[5] = 1.0;
        mat[6] = 0.0;
        mat[7] = 0.0;

        mat[8] = -sinAngle;
        mat[9] = 0.0;
        mat[10] = cosAngle;
        mat[11] = 0.0;

        mat[12] = 0.0;
        mat[13] = 0.0;
        mat[14] = 0.0;
        mat[15] = 1.0;

        // Issue 253: set all dirty bits if input is infinity or NaN
        if (isInfOrNaN(angle)) {
            dirtyBits = ALL_DIRTY;
            return;
        }

	type = CONGRUENT | AFFINE | RIGID | ORTHO;
	dirtyBits = CLASSIFY_BIT | ROTATION_BIT | SCALE_BIT;
    }


    /**
     * Sets the value of this transform to a counter clockwise rotation
     * about the z axis.  All of the non-rotational components are set
     * as if this were an identity matrix.
     * @param angle the angle to rotate about the Z axis in radians
     */
    public void rotZ(double angle)  {
        double sinAngle = Math.sin(angle);
        double cosAngle = Math.cos(angle);

        mat[0] = cosAngle;
        mat[1] = -sinAngle;
        mat[2] = 0.0;
        mat[3] = 0.0;

        mat[4] = sinAngle;
        mat[5] = cosAngle;
        mat[6] = 0.0;
        mat[7] = 0.0;

        mat[8] = 0.0;
        mat[9] = 0.0;
        mat[10] = 1.0;
        mat[11] = 0.0;

        mat[12] = 0.0;
        mat[13] = 0.0;
        mat[14] = 0.0;
        mat[15] = 1.0;

        // Issue 253: set all dirty bits if input is infinity or NaN
        if (isInfOrNaN(angle)) {
            dirtyBits = ALL_DIRTY;
            return;
        }

	type = CONGRUENT | AFFINE | RIGID | ORTHO;
	dirtyBits = CLASSIFY_BIT | ROTATION_BIT | SCALE_BIT;
    }

    /**
     * Sets the scale component of the current transform; any existing
     * scale is first factored out of the existing transform before
     * the new scale is applied.
     * @param scale  the new scale amount
     */
    public final void setScale(double scale) {
	if ((dirtyBits & ROTATION_BIT)!= 0) {
	    computeScaleRotation(false);
	}

	scales[0] = scales[1] = scales[2] = scale;
	mat[0] = rot[0]*scale;
	mat[1] = rot[1]*scale;
	mat[2] = rot[2]*scale;
	mat[4] = rot[3]*scale;
	mat[5] = rot[4]*scale;
	mat[6] = rot[5]*scale;
	mat[8] = rot[6]*scale;
	mat[9] = rot[7]*scale;
	mat[10] = rot[8]*scale;

        // Issue 253: set all dirty bits if input is infinity or NaN
        if (isInfOrNaN(scale)) {
            dirtyBits = ALL_DIRTY;
            return;
        }

	dirtyBits |= (CLASSIFY_BIT | RIGID_BIT | CONGRUENT_BIT | SVD_BIT);
	dirtyBits &= ~SCALE_BIT;
    }


    /**
     * Sets the possibly non-uniform scale component of the current
     * transform; any existing scale is first factored out of the
     * existing transform before the new scale is applied.
     * @param scale  the new x,y,z scale values
     */
     public final void setScale(IVecteur3d scale) {

	if ((dirtyBits & ROTATION_BIT)!= 0) {
	    computeScaleRotation(false);
	}

	scales[0] = scale.obtenirX();
	scales[1] = scale.obtenirY();
	scales[2] = scale.obtenirZ();

	mat[0] = rot[0]*scale.obtenirX();
	mat[1] = rot[1]*scale.obtenirY();
	mat[2] = rot[2]*scale.obtenirZ();
	mat[4] = rot[3]*scale.obtenirX();
	mat[5] = rot[4]*scale.obtenirY();
	mat[6] = rot[5]*scale.obtenirZ();
	mat[8] = rot[6]*scale.obtenirX();
	mat[9] = rot[7]*scale.obtenirY();
	mat[10] = rot[8]*scale.obtenirZ();

        // Issue 253: set all dirty bits if input is infinity or NaN
        if (isInfOrNaN(scale)) {
            dirtyBits = ALL_DIRTY;
            return;
        }

	dirtyBits |= (CLASSIFY_BIT | RIGID_BIT | CONGRUENT_BIT | SVD_BIT);
	dirtyBits &= ~SCALE_BIT;
    }


    /**
     * Replaces the current transform with a non-uniform scale transform.
     * All values of the existing transform are replaced.
     * @param xScale the new X scale amount
     * @param yScale the new Y scale amount
     * @param zScale the new Z scale amount
     * @deprecated Use setScale(Vector3d) instead of setNonUniformScale;
     * note that the setScale only modifies the scale component
     */
    @Deprecated
	public final void setNonUniformScale(double xScale,
					 double yScale,
					 double zScale) {
	if(scales == null)
	    scales = new double[3];

	scales[0] = xScale;
	scales[1] = yScale;
	scales[2] = zScale;
	mat[0] = xScale;
	mat[1] = 0.0;
	mat[2] = 0.0;
	mat[3] = 0.0;
	mat[4] = 0.0;
	mat[5] = yScale;
	mat[6] = 0.0;
	mat[7] = 0.0;
	mat[8] = 0.0;
	mat[9] = 0.0;
	mat[10] = zScale;
	mat[11] = 0.0;
	mat[12] = 0.0;
	mat[13] = 0.0;
	mat[14] = 0.0;
	mat[15] = 1.0;

        // Issue 253: set all dirty bits
        dirtyBits = ALL_DIRTY;
    }

    /**
     * Replaces the translational components of this transform to the values
     * in the Vector3f argument; the other values of this transform are not
     * modified.
     * @param trans  the translational component
     */
    public final void setTranslation(IVecteur3d trans) {
       mat[3] = trans.obtenirX();
       mat[7] = trans.obtenirY();
       mat[11] = trans.obtenirZ();

       // Issue 253: set all dirty bits if input is infinity or NaN
       if (isInfOrNaN(trans)) {
           dirtyBits = ALL_DIRTY;
           return;
       }

       // Only preserve CONGRUENT, RIGID, ORTHO
       type &= ~(ORTHOGONAL|IDENTITY|SCALE|TRANSLATION|SCALE|ZERO);
       dirtyBits |= CLASSIFY_BIT;
    }


    /**
     * Sets the rotational component (upper 3x3) of this transform to the
     * rotation matrix converted from the Euler angles provided; the other
     * non-rotational elements are set as if this were an identity matrix.
     * The euler parameter is a Vector3d consisting of three rotation angles
     * applied first about the X, then Y then Z axis.
     * These rotations are applied using a static frame of reference. In
     * other words, the orientation of the Y rotation axis is not affected
     * by the X rotation and the orientation of the Z rotation axis is not
     * affected by the X or Y rotation.
     * @param euler  the Vector3d consisting of three rotation angles about X,Y,Z
     *
     */
    public final void setEuler(IVecteur3d euler) {
	double sina, sinb, sinc;
	double cosa, cosb, cosc;

	sina = Math.sin(euler.obtenirX());
	sinb = Math.sin(euler.obtenirY());
	sinc = Math.sin(euler.obtenirZ());
	cosa = Math.cos(euler.obtenirX());
	cosb = Math.cos(euler.obtenirY());
	cosc = Math.cos(euler.obtenirZ());

	mat[0] = cosb * cosc;
	mat[1] = -(cosa * sinc) + (sina * sinb * cosc);
	mat[2] = (sina * sinc) + (cosa * sinb *cosc);
	mat[3] = 0.0;

	mat[4] = cosb * sinc;
	mat[5] = (cosa * cosc) + (sina * sinb * sinc);
	mat[6] = -(sina * cosc) + (cosa * sinb *sinc);
	mat[7] = 0.0;

	mat[8] = -sinb;
	mat[9] = sina * cosb;
	mat[10] = cosa * cosb;
	mat[11] = 0.0;

	mat[12] = 0.0;
	mat[13] = 0.0;
	mat[14] = 0.0;
	mat[15] = 1.0;

        // Issue 253: set all dirty bits if input is infinity or NaN
        if (isInfOrNaN(euler)) {
            dirtyBits = ALL_DIRTY;
            return;
        }

	type = AFFINE | CONGRUENT | RIGID | ORTHO;
	dirtyBits = CLASSIFY_BIT | SCALE_BIT | ROTATION_BIT;
    }
    /**
     * Returns the uniform scale factor of this matrix.
     * If the matrix has non-uniform scale factors, the largest of the
     * x, y, and z scale factors will be returned.
     * @return  the scale factor of this matrix
     */
    public final double getScale() {
	if ((dirtyBits & SCALE_BIT) != 0) {
	    computeScales(false);
	}
	return max3(scales);
   }


    /**
     * Gets the possibly non-uniform scale components of the current
     * transform and places them into the scale vector.
     * @param scale  the vector into which the x,y,z scale values will be placed
     */
    public final void getScale(IVecteur3d scale) {
	if ((dirtyBits & SCALE_BIT) != 0) {
	    computeScales(false);
	}
	scale.définirX(scales[0]);
	scale.définirY(scales[1]);
	scale.définirZ(scales[2]);
    }

    /**
     * Sets the value of this transform to the inverse of the passed
     * Transformation3D parameter.  This method uses the transform type
     * to determine the optimal algorithm for inverting transform t1.
     * @param t1  the transform to be inverted
     * @exception SingularMatrixException thrown if transform t1 is
     * not invertible
     */
    public final void invert(Transformation3D t1) {
	if (t1 == this) {
	    invert();
	} else if (t1.isAffine()) {
	    // We can't use invertOrtho() because of numerical
	    // instability unless we set tolerance of ortho test to 0
	    invertAffine(t1);
	} else {
	    invertGeneral(t1);
	}
    }

    /**
     * Inverts this transform in place.  This method uses the transform
     * type to determine the optimal algorithm for inverting this transform.
     * @exception SingularMatrixException thrown if this transform is
     * not invertible
     */
    public final void invert() {
	if (isAffine()) {
	    invertAffine();
	} else {
	    invertGeneral(this);
	}
    }

    /**
     * Affine invert routine.  Inverts t1 and places the result in "this".
     */
    final void invertAffine(Transformation3D t1) {
	double determinant = t1.affineDeterminant();


	double s = (t1.mat[0]*t1.mat[0] + t1.mat[1]*t1.mat[1] +
		    t1.mat[2]*t1.mat[2] + t1.mat[3]*t1.mat[3])*
	           (t1.mat[4]*t1.mat[4] + t1.mat[5]*t1.mat[5] +
	            t1.mat[6]*t1.mat[6] + t1.mat[7]*t1.mat[7])*
	           (t1.mat[8]*t1.mat[8] + t1.mat[9]*t1.mat[9] +
		    t1.mat[10]*t1.mat[10] + t1.mat[11]*t1.mat[11]);

	if ((determinant*determinant) < (EPS * s)) {
	    // using invertGeneral is numerically more stable for
	    //this case  see bug 4227733
	    invertGeneral(t1);
	    return;
	}
	s = 1.0 / determinant;

	mat[0] =  (t1.mat[5]*t1.mat[10] - t1.mat[9]*t1.mat[6]) * s;
	mat[1] = -(t1.mat[1]*t1.mat[10] - t1.mat[9]*t1.mat[2]) * s;
	mat[2] =  (t1.mat[1]*t1.mat[6] - t1.mat[5]*t1.mat[2]) * s;
	mat[4] = -(t1.mat[4]*t1.mat[10] - t1.mat[8]*t1.mat[6]) * s;
	mat[5] =  (t1.mat[0]*t1.mat[10] - t1.mat[8]*t1.mat[2]) * s;
	mat[6] = -(t1.mat[0]*t1.mat[6] -  t1.mat[4]*t1.mat[2]) * s;
	mat[8] =  (t1.mat[4]*t1.mat[9] - t1.mat[8]*t1.mat[5]) * s;
	mat[9] = -(t1.mat[0]*t1.mat[9] - t1.mat[8]*t1.mat[1]) * s;
	mat[10]=  (t1.mat[0]*t1.mat[5] - t1.mat[4]*t1.mat[1]) * s;
	mat[3] = -(t1.mat[3] * mat[0] + t1.mat[7] * mat[1] +
		   t1.mat[11] * mat[2]);
	mat[7] = -(t1.mat[3] * mat[4] + t1.mat[7] * mat[5] +
		   t1.mat[11] * mat[6]);
	mat[11]= -(t1.mat[3] * mat[8] + t1.mat[7] * mat[9] +
		   t1.mat[11] * mat[10]);

	mat[12] = mat[13] = mat[14] = 0.0;
	mat[15] = 1.0;

	dirtyBits = t1.dirtyBits | ROTSCALESVD_DIRTY | CLASSIFY_BIT | ORTHO_BIT;
	type = t1.type;
    }

    /**
     * Affine invert routine.  Inverts "this" matrix in place.
     */
    final void invertAffine() {
	double determinant = affineDeterminant();

	double s = (mat[0]*mat[0] + mat[1]*mat[1] +
	            mat[2]*mat[2] + mat[3]*mat[3])*
	           (mat[4]*mat[4] + mat[5]*mat[5] +
	            mat[6]*mat[6] + mat[7]*mat[7])*
	           (mat[8]*mat[8] + mat[9]*mat[9] +
		    mat[10]*mat[10] + mat[11]*mat[11]);

	if ((determinant*determinant) < (EPS * s)) {
	    invertGeneral(this);
	    return;
	}
	s = 1.0 / determinant;
	double tmp0 =  (mat[5]*mat[10] - mat[9]*mat[6]) * s;
	double tmp1 = -(mat[1]*mat[10] - mat[9]*mat[2]) * s;
	double tmp2 =  (mat[1]*mat[6] -  mat[5]*mat[2]) * s;
	double tmp4 = -(mat[4]*mat[10] - mat[8]*mat[6]) * s;
	double tmp5 =  (mat[0]*mat[10] - mat[8]*mat[2]) * s;
	double tmp6 = -(mat[0]*mat[6]  - mat[4]*mat[2]) * s;
	double tmp8 =  (mat[4]*mat[9]  - mat[8]*mat[5]) * s;
	double tmp9 = -(mat[0]*mat[9]  - mat[8]*mat[1]) * s;
	double tmp10=  (mat[0]*mat[5]  - mat[4]*mat[1]) * s;
	double tmp3 = -(mat[3] * tmp0 + mat[7] * tmp1 + mat[11] * tmp2);
	double tmp7 = -(mat[3] * tmp4 + mat[7] * tmp5 + mat[11] * tmp6);
	mat[11]= -(mat[3] * tmp8 + mat[7] * tmp9 + mat[11] * tmp10);

	mat[0]=tmp0; mat[1]=tmp1; mat[2]=tmp2; mat[3]=tmp3;
	mat[4]=tmp4; mat[5]=tmp5; mat[6]=tmp6; mat[7]=tmp7;
	mat[8]=tmp8; mat[9]=tmp9; mat[10]=tmp10;
	mat[12] = mat[13] = mat[14] = 0.0;
	mat[15] = 1.0;
	dirtyBits |= ROTSCALESVD_DIRTY | CLASSIFY_BIT | ORTHO_BIT;
    }

    /**
     * General invert routine.  Inverts t1 and places the result in "this".
     * Note that this routine handles both the "this" version and the
     * non-"this" version.
     *
     * Also note that since this routine is slow anyway, we won't worry
     * about allocating a little bit of garbage.
     */
    final void invertGeneral(Transformation3D t1) {
	double tmp[] = new double[16];
	int row_perm[] = new int[4];

	// Use LU decomposition and backsubstitution code specifically
	// for floating-point 4x4 matrices.

	// Copy source matrix to tmp
	System.arraycopy(t1.mat, 0, tmp, 0, tmp.length);

	// Perform back substitution on the identity matrix
	// luDecomposition will set rot[] & scales[] for use
	// in luBacksubstituation
	mat[0] = 1.0;  mat[1] = 0.0;  mat[2] = 0.0;  mat[3] = 0.0;
	mat[4] = 0.0;  mat[5] = 1.0;  mat[6] = 0.0;  mat[7] = 0.0;
	mat[8] = 0.0;  mat[9] = 0.0;  mat[10] = 1.0; mat[11] = 0.0;
	mat[12] = 0.0; mat[13] = 0.0; mat[14] = 0.0; mat[15] = 1.0;
	luBacksubstitution(tmp, row_perm, this.mat);

	type = 0;
	dirtyBits = ALL_DIRTY;
    }


    /**
     * Given a 4x4 array "matrix0", this function replaces it with the
     * LU decomposition of a row-wise permutation of itself.  The input
     * parameters are "matrix0" and "dimen".  The array "matrix0" is also
     * an output parameter.  The vector "row_perm[4]" is an output
     * parameter that contains the row permutations resulting from partial
     * pivoting.  The output parameter "even_row_xchg" is 1 when the
     * number of row exchanges is even, or -1 otherwise.  Assumes data
     * type is always double.
     *
     * This function is similar to luDecomposition, except that it
     * is tuned specifically for 4x4 matrices.
     *
     * @return true if the matrix is nonsingular, or false otherwise.
     */
    //
    // Reference: Press, Flannery, Teukolsky, Vetterling,
    //	      _Numerical_Recipes_in_C_, Cambridge University Press,
    //	      1988, pp 40-45.
    //
    static boolean luDecomposition(double[] matrix0,
				   int[] row_perm) {

	// Can't re-use this temporary since the method is static.
	double row_scale[] = new double[4];

	// Determine implicit scaling information by looping over rows
	{
	    int i, j;
	    int ptr, rs;
	    double big, temp;

	    ptr = 0;
	    rs = 0;

	    // For each row ...
	    i = 4;
	    while (i-- != 0) {
		big = 0.0;

		// For each column, find the largest element in the row
		j = 4;
		while (j-- != 0) {
		    temp = matrix0[ptr++];
		    temp = Math.abs(temp);
		    if (temp > big) {
			big = temp;
		    }
		}

		// Is the matrix singular?
		if (big == 0.0) {
		    return false;
		}
		row_scale[rs++] = 1.0 / big;
	    }
	}

	{
	    int j;
	    int mtx;

	    mtx = 0;

	    // For all columns, execute Crout's method
	    for (j = 0; j < 4; j++) {
		int i, imax, k;
		int target, p1, p2;
		double sum, big, temp;

		// Determine elements of upper diagonal matrix U
		for (i = 0; i < j; i++) {
		    target = mtx + (4*i) + j;
		    sum = matrix0[target];
		    k = i;
		    p1 = mtx + (4*i);
		    p2 = mtx + j;
		    while (k-- != 0) {
			sum -= matrix0[p1] * matrix0[p2];
			p1++;
			p2 += 4;
		    }
		    matrix0[target] = sum;
		}

		// Search for largest pivot element and calculate
		// intermediate elements of lower diagonal matrix L.
		big = 0.0;
		imax = -1;
		for (i = j; i < 4; i++) {
		    target = mtx + (4*i) + j;
		    sum = matrix0[target];
		    k = j;
		    p1 = mtx + (4*i);
		    p2 = mtx + j;
		    while (k-- != 0) {
			sum -= matrix0[p1] * matrix0[p2];
			p1++;
			p2 += 4;
		    }
		    matrix0[target] = sum;

		    // Is this the best pivot so far?
		    if ((temp = row_scale[i] * Math.abs(sum)) >= big) {
			big = temp;
			imax = i;
		    }
		}

		if (imax < 0) {
		    return false;
		}

		// Is a row exchange necessary?
		if (j != imax) {
		    // Yes: exchange rows
		    k = 4;
		    p1 = mtx + (4*imax);
		    p2 = mtx + (4*j);
		    while (k-- != 0) {
			temp = matrix0[p1];
			matrix0[p1++] = matrix0[p2];
			matrix0[p2++] = temp;
		    }

		    // Record change in scale factor
		    row_scale[imax] = row_scale[j];
		}

		// Record row permutation
		row_perm[j] = imax;

		// Is the matrix singular
		if (matrix0[(mtx + (4*j) + j)] == 0.0) {
		    return false;
		}

		// Divide elements of lower diagonal matrix L by pivot
		if (j != (4-1)) {
		    temp = 1.0 / (matrix0[(mtx + (4*j) + j)]);
		    target = mtx + (4*(j+1)) + j;
		    i = 3 - j;
		    while (i-- != 0) {
			matrix0[target] *= temp;
			target += 4;
		    }
		}
	    }
	}

	return true;
    }


    /**
     * Solves a set of linear equations.  The input parameters "matrix1",
     * and "row_perm" come from luDecompostionD4x4 and do not change
     * here.  The parameter "matrix2" is a set of column vectors assembled
     * into a 4x4 matrix of floating-point values.  The procedure takes each
     * column of "matrix2" in turn and treats it as the right-hand side of the
     * matrix equation Ax = LUx = b.  The solution vector replaces the
     * original column of the matrix.
     *
     * If "matrix2" is the identity matrix, the procedure replaces its contents
     * with the inverse of the matrix from which "matrix1" was originally
     * derived.
     */
    //
    // Reference: Press, Flannery, Teukolsky, Vetterling,
    //	      _Numerical_Recipes_in_C_, Cambridge University Press,
    //	      1988, pp 44-45.
    //
    static void luBacksubstitution(double[] matrix1,
				   int[] row_perm,
				   double[] matrix2) {

	int i, ii, ip, j, k;
	int rp;
	int cv, rv;

	//	rp = row_perm;
	rp = 0;

	// For each column vector of matrix2 ...
	for (k = 0; k < 4; k++) {
	    //	    cv = &(matrix2[0][k]);
	    cv = k;
	    ii = -1;

	    // Forward substitution
	    for (i = 0; i < 4; i++) {
		double sum;

		ip = row_perm[rp+i];
		sum = matrix2[cv+4*ip];
		matrix2[cv+4*ip] = matrix2[cv+4*i];
		if (ii >= 0) {
		    //		    rv = &(matrix1[i][0]);
		    rv = i*4;
		    for (j = ii; j <= i-1; j++) {
			sum -= matrix1[rv+j] * matrix2[cv+4*j];
		    }
		}
		else if (sum != 0.0) {
		    ii = i;
		}
		matrix2[cv+4*i] = sum;
	    }

	    // Backsubstitution
	    //	    rv = &(matrix1[3][0]);
	    rv = 3*4;
	    matrix2[cv+4*3] /= matrix1[rv+3];

	    rv -= 4;
	    matrix2[cv+4*2] = (matrix2[cv+4*2] -
			    matrix1[rv+3] * matrix2[cv+4*3]) / matrix1[rv+2];

	    rv -= 4;
	    matrix2[cv+4*1] = (matrix2[cv+4*1] -
			    matrix1[rv+2] * matrix2[cv+4*2] -
			    matrix1[rv+3] * matrix2[cv+4*3]) / matrix1[rv+1];

	    rv -= 4;
	    matrix2[cv+4*0] = (matrix2[cv+4*0] -
			    matrix1[rv+1] * matrix2[cv+4*1] -
			    matrix1[rv+2] * matrix2[cv+4*2] -
			    matrix1[rv+3] * matrix2[cv+4*3]) / matrix1[rv+0];
	}
    }

    // given that this matrix is affine
    final double affineDeterminant() {
	return mat[0]*(mat[5]*mat[10] - mat[6]*mat[9]) -
	       mat[1]*(mat[4]*mat[10] - mat[6]*mat[8]) +
	       mat[2]*(mat[4]*mat[ 9] - mat[5]*mat[8]);
    }

    /**
     * Calculates and returns the determinant of this transform.
     * @return  the double precision determinant
     */
     public final double determinant() {

	 if (isAffine()) {
	     return mat[0]*(mat[5]*mat[10] - mat[6]*mat[9]) -
	 	    mat[1]*(mat[4]*mat[10] - mat[6]*mat[8]) +
		    mat[2]*(mat[4]*mat[ 9] - mat[5]*mat[8]);
	 }
	 // cofactor exapainsion along first row
	 return mat[0]*(mat[5]*(mat[10]*mat[15] - mat[11]*mat[14]) -
			mat[6]*(mat[ 9]*mat[15] - mat[11]*mat[13]) +
			mat[7]*(mat[ 9]*mat[14] - mat[10]*mat[13])) -
	        mat[1]*(mat[4]*(mat[10]*mat[15] - mat[11]*mat[14]) -
	                mat[6]*(mat[ 8]*mat[15] - mat[11]*mat[12]) +
			mat[7]*(mat[ 8]*mat[14] - mat[10]*mat[12])) +
	        mat[2]*(mat[4]*(mat[ 9]*mat[15] - mat[11]*mat[13]) -
			mat[5]*(mat[ 8]*mat[15] - mat[11]*mat[12]) +
			mat[7]*(mat[ 8]*mat[13] - mat[ 9]*mat[12])) -
	        mat[3]*(mat[4]*(mat[ 9]*mat[14] - mat[10]*mat[13]) -
			mat[5]*(mat[ 8]*mat[14] - mat[10]*mat[12]) +
			mat[6]*(mat[ 8]*mat[13] - mat[ 9]*mat[12]));
     }



   /**
     * Multiplies each element of this transform by a scalar.
     * @param scalar  the scalar multiplier
     */
    public final void mul(double scalar) {
	for (int i=0 ; i<16 ; i++) {
	    mat[i] *= scalar;
	}
	dirtyBits = ALL_DIRTY;
    }

   /**
     * Multiplies each element of transform t1 by a scalar and places
     * the result into this.  Transform t1 is not modified.
     * @param scalar  the scalar multiplier
     * @param t1  the original transform
     */
    public final void mul(double scalar, Transformation3D t1)  {
	for (int i=0 ; i<16 ; i++) {
	    mat[i] = t1.mat[i] * scalar;
	}
	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}

    }


    /**
     * Sets the value of this transform to the result of multiplying itself
     * with transform t1 (this = this * t1).
     * @param t1 the other transform
     */
    public final void mul(Transformation3D t1) {
	double tmp0, tmp1, tmp2, tmp3;
	double tmp4, tmp5, tmp6, tmp7;
	double tmp8, tmp9, tmp10, tmp11;
	boolean aff = false;

	if (t1.isAffine()) {
	    tmp0 = mat[0]*t1.mat[0] + mat[1]*t1.mat[4] + mat[2]*t1.mat[8];
	    tmp1 = mat[0]*t1.mat[1] + mat[1]*t1.mat[5] + mat[2]*t1.mat[9];
	    tmp2 = mat[0]*t1.mat[2] + mat[1]*t1.mat[6] + mat[2]*t1.mat[10];
	    tmp3 = mat[0]*t1.mat[3] + mat[1]*t1.mat[7] + mat[2]*t1.mat[11] + mat[3];
	    tmp4 = mat[4]*t1.mat[0] + mat[5]*t1.mat[4] + mat[6]*t1.mat[8];
	    tmp5 = mat[4]*t1.mat[1] + mat[5]*t1.mat[5] + mat[6]*t1.mat[9];
	    tmp6 = mat[4]*t1.mat[2] + mat[5]*t1.mat[6] + mat[6]*t1.mat[10];
	    tmp7 = mat[4]*t1.mat[3] + mat[5]*t1.mat[7] + mat[6]*t1.mat[11] + mat[7];
	    tmp8 = mat[8]*t1.mat[0] + mat[9]*t1.mat[4] + mat[10]*t1.mat[8];
	    tmp9 = mat[8]*t1.mat[1] + mat[9]*t1.mat[5] + mat[10]*t1.mat[9];
	    tmp10 = mat[8]*t1.mat[2] + mat[9]*t1.mat[6] + mat[10]*t1.mat[10];
	    tmp11 = mat[8]*t1.mat[3] + mat[9]*t1.mat[7] + mat[10]*t1.mat[11] + mat[11];
	    if (isAffine()) {
		mat[12] =  mat[13] = mat[14] = 0;
		mat[15] = 1;
		aff = true;
	    } else {
		double tmp12 = mat[12]*t1.mat[0] + mat[13]*t1.mat[4] +
		               mat[14]*t1.mat[8];
		double tmp13 = mat[12]*t1.mat[1] + mat[13]*t1.mat[5] +
		               mat[14]*t1.mat[9];
		double tmp14 = mat[12]*t1.mat[2] + mat[13]*t1.mat[6] +
		               mat[14]*t1.mat[10];
                double tmp15 = mat[12]*t1.mat[3] + mat[13]*t1.mat[7] +
		               mat[14]*t1.mat[11] + mat[15];
		mat[12] = tmp12;
		mat[13] = tmp13;
		mat[14] = tmp14;
		mat[15] = tmp15;
	    }
	} else {
	    tmp0 = mat[0]*t1.mat[0] + mat[1]*t1.mat[4] + mat[2]*t1.mat[8] +
                   mat[3]*t1.mat[12];
	    tmp1 = mat[0]*t1.mat[1] + mat[1]*t1.mat[5] + mat[2]*t1.mat[9] +
                   mat[3]*t1.mat[13];
	    tmp2 = mat[0]*t1.mat[2] + mat[1]*t1.mat[6] + mat[2]*t1.mat[10] +
                   mat[3]*t1.mat[14];
	    tmp3 = mat[0]*t1.mat[3] + mat[1]*t1.mat[7] + mat[2]*t1.mat[11] +
                   mat[3]*t1.mat[15];
            tmp4 = mat[4]*t1.mat[0] + mat[5]*t1.mat[4] + mat[6]*t1.mat[8] +
		   mat[7]*t1.mat[12];
	    tmp5 = mat[4]*t1.mat[1] + mat[5]*t1.mat[5] + mat[6]*t1.mat[9] +
                   mat[7]*t1.mat[13];
	    tmp6 = mat[4]*t1.mat[2] + mat[5]*t1.mat[6] + mat[6]*t1.mat[10] +
                   mat[7]*t1.mat[14];
	    tmp7 = mat[4]*t1.mat[3] + mat[5]*t1.mat[7] + mat[6]*t1.mat[11] +
                   mat[7]*t1.mat[15];
            tmp8 = mat[8]*t1.mat[0] + mat[9]*t1.mat[4] + mat[10]*t1.mat[8] +
                   mat[11]*t1.mat[12];
            tmp9 = mat[8]*t1.mat[1] + mat[9]*t1.mat[5] + mat[10]*t1.mat[9] +
                   mat[11]*t1.mat[13];
            tmp10 = mat[8]*t1.mat[2] + mat[9]*t1.mat[6] +
		    mat[10]*t1.mat[10]+ mat[11]*t1.mat[14];
            tmp11 = mat[8]*t1.mat[3] + mat[9]*t1.mat[7] +
		    mat[10]*t1.mat[11] + mat[11]*t1.mat[15];

	    if (isAffine()) {
		mat[12] = t1.mat[12];
		mat[13] = t1.mat[13];
		mat[14] = t1.mat[14];
		mat[15] = t1.mat[15];
	    } else {
		double tmp12 = mat[12]*t1.mat[0] + mat[13]*t1.mat[4] +
		               mat[14]*t1.mat[8] +  mat[15]*t1.mat[12];
		double tmp13 = mat[12]*t1.mat[1] + mat[13]*t1.mat[5] +
		               mat[14]*t1.mat[9] + mat[15]*t1.mat[13];
		double tmp14 = mat[12]*t1.mat[2] + mat[13]*t1.mat[6] +
		               mat[14]*t1.mat[10] + mat[15]*t1.mat[14];
		double tmp15 = mat[12]*t1.mat[3] + mat[13]*t1.mat[7] +
		               mat[14]*t1.mat[11] + mat[15]*t1.mat[15];
		mat[12] = tmp12;
		mat[13] = tmp13;
		mat[14] = tmp14;
		mat[15] = tmp15;
	    }
	}

	mat[0] = tmp0;
	mat[1] = tmp1;
	mat[2] = tmp2;
	mat[3] = tmp3;
	mat[4] = tmp4;
	mat[5] = tmp5;
	mat[6] = tmp6;
	mat[7] = tmp7;
	mat[8] = tmp8;
	mat[9] = tmp9;
	mat[10] = tmp10;
	mat[11] = tmp11;

	if (((dirtyBits & CONGRUENT_BIT) == 0) &&
	    ((type & CONGRUENT) != 0) &&
	    ((t1.dirtyBits & CONGRUENT_BIT) == 0) &&
	    ((t1.type & CONGRUENT) != 0)) {
	    type &= t1.type;
	    dirtyBits |= t1.dirtyBits | CLASSIFY_BIT |
		ROTSCALESVD_DIRTY | RIGID_BIT;
	} else {
	    if (aff) {
		dirtyBits = ORTHO_BIT | CONGRUENT_BIT | RIGID_BIT |
                            CLASSIFY_BIT | ROTSCALESVD_DIRTY;
	    } else {
		dirtyBits = ALL_DIRTY;
	    }
	}

	if (autoNormalize) {
	    normalize();
	}

    }

    /**
     * Sets the value of this transform to the result of multiplying transform
     * t1 by transform t2 (this = t1*t2).
     * @param t1  the left transform
     * @param t2  the right transform
     */
    public final void mul(Transformation3D t1, Transformation3D t2) {
	boolean aff = false;
	if ((this != t1)  &&  (this != t2)) {
	    if (t2.isAffine()) {

		mat[0] = t1.mat[0]*t2.mat[0] + t1.mat[1]*t2.mat[4] + t1.mat[2]*t2.mat[8];
		mat[1] = t1.mat[0]*t2.mat[1] + t1.mat[1]*t2.mat[5] + t1.mat[2]*t2.mat[9];
		mat[2] = t1.mat[0]*t2.mat[2] + t1.mat[1]*t2.mat[6] + t1.mat[2]*t2.mat[10];
		mat[3] = t1.mat[0]*t2.mat[3] + t1.mat[1]*t2.mat[7] +
                         t1.mat[2]*t2.mat[11] + t1.mat[3];
		mat[4] = t1.mat[4]*t2.mat[0] + t1.mat[5]*t2.mat[4] + t1.mat[6]*t2.mat[8];
		mat[5] = t1.mat[4]*t2.mat[1] + t1.mat[5]*t2.mat[5] + t1.mat[6]*t2.mat[9];
	        mat[6] = t1.mat[4]*t2.mat[2] + t1.mat[5]*t2.mat[6] + t1.mat[6]*t2.mat[10];
	        mat[7] = t1.mat[4]*t2.mat[3] + t1.mat[5]*t2.mat[7] +
                         t1.mat[6]*t2.mat[11] + t1.mat[7];
		mat[8] = t1.mat[8]*t2.mat[0] + t1.mat[9]*t2.mat[4] + t1.mat[10]*t2.mat[8];
	        mat[9] = t1.mat[8]*t2.mat[1] + t1.mat[9]*t2.mat[5] + t1.mat[10]*t2.mat[9];
		mat[10] = t1.mat[8]*t2.mat[2] + t1.mat[9]*t2.mat[6] + t1.mat[10]*t2.mat[10];
	        mat[11] = t1.mat[8]*t2.mat[3] + t1.mat[9]*t2.mat[7] +
                          t1.mat[10]*t2.mat[11] + t1.mat[11];
		if (t1.isAffine()) {
		    aff = true;
		    mat[12] =  mat[13] = mat[14] = 0;
		    mat[15] = 1;
		} else {
		    mat[12] = t1.mat[12]*t2.mat[0] + t1.mat[13]*t2.mat[4] +
			      t1.mat[14]*t2.mat[8];
		    mat[13] = t1.mat[12]*t2.mat[1] + t1.mat[13]*t2.mat[5] +
		              t1.mat[14]*t2.mat[9];
		    mat[14] = t1.mat[12]*t2.mat[2] + t1.mat[13]*t2.mat[6] +
		              t1.mat[14]*t2.mat[10];
		    mat[15] = t1.mat[12]*t2.mat[3] + t1.mat[13]*t2.mat[7] +
			      t1.mat[14]*t2.mat[11] + t1.mat[15];
		}
	    } else {
		mat[0] = t1.mat[0]*t2.mat[0] + t1.mat[1]*t2.mat[4] +
		         t1.mat[2]*t2.mat[8] + t1.mat[3]*t2.mat[12];
		mat[1] = t1.mat[0]*t2.mat[1] + t1.mat[1]*t2.mat[5] +
		         t1.mat[2]*t2.mat[9] + t1.mat[3]*t2.mat[13];
		mat[2] = t1.mat[0]*t2.mat[2] + t1.mat[1]*t2.mat[6] +
		         t1.mat[2]*t2.mat[10] + t1.mat[3]*t2.mat[14];
		mat[3] = t1.mat[0]*t2.mat[3] + t1.mat[1]*t2.mat[7] +
		         t1.mat[2]*t2.mat[11] + t1.mat[3]*t2.mat[15];
		mat[4] = t1.mat[4]*t2.mat[0] + t1.mat[5]*t2.mat[4] +
		         t1.mat[6]*t2.mat[8] + t1.mat[7]*t2.mat[12];
		mat[5] = t1.mat[4]*t2.mat[1] + t1.mat[5]*t2.mat[5] +
		         t1.mat[6]*t2.mat[9] + t1.mat[7]*t2.mat[13];
		mat[6] = t1.mat[4]*t2.mat[2] + t1.mat[5]*t2.mat[6] +
		         t1.mat[6]*t2.mat[10] + t1.mat[7]*t2.mat[14];
		mat[7] = t1.mat[4]*t2.mat[3] + t1.mat[5]*t2.mat[7] +
		         t1.mat[6]*t2.mat[11] + t1.mat[7]*t2.mat[15];
		mat[8] = t1.mat[8]*t2.mat[0] + t1.mat[9]*t2.mat[4] +
		         t1.mat[10]*t2.mat[8] + t1.mat[11]*t2.mat[12];
		mat[9] = t1.mat[8]*t2.mat[1] + t1.mat[9]*t2.mat[5] +
		         t1.mat[10]*t2.mat[9] + t1.mat[11]*t2.mat[13];
		mat[10] = t1.mat[8]*t2.mat[2] + t1.mat[9]*t2.mat[6] +
		          t1.mat[10]*t2.mat[10] + t1.mat[11]*t2.mat[14];
		mat[11] = t1.mat[8]*t2.mat[3] + t1.mat[9]*t2.mat[7] +
		          t1.mat[10]*t2.mat[11] + t1.mat[11]*t2.mat[15];
		if (t1.isAffine()) {
		    mat[12] = t2.mat[12];
		    mat[13] = t2.mat[13];
		    mat[14] = t2.mat[14];
		    mat[15] = t2.mat[15];
		} else {
		    mat[12] = t1.mat[12]*t2.mat[0] + t1.mat[13]*t2.mat[4] +
		              t1.mat[14]*t2.mat[8] + t1.mat[15]*t2.mat[12];
		    mat[13] = t1.mat[12]*t2.mat[1] + t1.mat[13]*t2.mat[5] +
		              t1.mat[14]*t2.mat[9] + t1.mat[15]*t2.mat[13];
		    mat[14] = t1.mat[12]*t2.mat[2] + t1.mat[13]*t2.mat[6] +
			      t1.mat[14]*t2.mat[10] + t1.mat[15]*t2.mat[14];
		    mat[15] = t1.mat[12]*t2.mat[3] + t1.mat[13]*t2.mat[7] +
		              t1.mat[14]*t2.mat[11] + t1.mat[15]*t2.mat[15];
		}
	    }
	} else {
	    double tmp0, tmp1, tmp2, tmp3;
	    double tmp4, tmp5, tmp6, tmp7;
	    double tmp8, tmp9, tmp10, tmp11;

	    if (t2.isAffine()) {
		tmp0 = t1.mat[0]*t2.mat[0] + t1.mat[1]*t2.mat[4] + t1.mat[2]*t2.mat[8];
		tmp1 = t1.mat[0]*t2.mat[1] + t1.mat[1]*t2.mat[5] + t1.mat[2]*t2.mat[9];
		tmp2 = t1.mat[0]*t2.mat[2] + t1.mat[1]*t2.mat[6] + t1.mat[2]*t2.mat[10];
		tmp3 = t1.mat[0]*t2.mat[3] + t1.mat[1]*t2.mat[7] +
		       t1.mat[2]*t2.mat[11] + t1.mat[3];
		tmp4 = t1.mat[4]*t2.mat[0] + t1.mat[5]*t2.mat[4] + t1.mat[6]*t2.mat[8];
		tmp5 = t1.mat[4]*t2.mat[1] + t1.mat[5]*t2.mat[5] + t1.mat[6]*t2.mat[9];
	        tmp6 = t1.mat[4]*t2.mat[2] + t1.mat[5]*t2.mat[6] + t1.mat[6]*t2.mat[10];
	        tmp7 = t1.mat[4]*t2.mat[3] + t1.mat[5]*t2.mat[7] +
                       t1.mat[6]*t2.mat[11] + t1.mat[7];
		tmp8 = t1.mat[8]*t2.mat[0] + t1.mat[9]*t2.mat[4] + t1.mat[10]*t2.mat[8];
	        tmp9 = t1.mat[8]*t2.mat[1] + t1.mat[9]*t2.mat[5] + t1.mat[10]*t2.mat[9];
		tmp10 = t1.mat[8]*t2.mat[2] + t1.mat[9]*t2.mat[6] + t1.mat[10]*t2.mat[10];
	        tmp11 = t1.mat[8]*t2.mat[3] + t1.mat[9]*t2.mat[7] +
                        t1.mat[10]*t2.mat[11] + t1.mat[11];
		if (t1.isAffine()) {
		    aff = true;
		    mat[12] =  mat[13] = mat[14] = 0;
		    mat[15] = 1;
		} else {
		    double tmp12 = t1.mat[12]*t2.mat[0] + t1.mat[13]*t2.mat[4] +
			           t1.mat[14]*t2.mat[8];
		    double tmp13 = t1.mat[12]*t2.mat[1] + t1.mat[13]*t2.mat[5] +
		                   t1.mat[14]*t2.mat[9];
		    double tmp14 = t1.mat[12]*t2.mat[2] + t1.mat[13]*t2.mat[6] +
		                   t1.mat[14]*t2.mat[10];
		    double tmp15 = t1.mat[12]*t2.mat[3] + t1.mat[13]*t2.mat[7] +
			           t1.mat[14]*t2.mat[11] + t1.mat[15];
		    mat[12] = tmp12;
		    mat[13] = tmp13;
		    mat[14] = tmp14;
		    mat[15] = tmp15;
		}
	    } else {
		tmp0 = t1.mat[0]*t2.mat[0] + t1.mat[1]*t2.mat[4] +
		       t1.mat[2]*t2.mat[8] + t1.mat[3]*t2.mat[12];
		tmp1 = t1.mat[0]*t2.mat[1] + t1.mat[1]*t2.mat[5] +
		       t1.mat[2]*t2.mat[9] + t1.mat[3]*t2.mat[13];
		tmp2 = t1.mat[0]*t2.mat[2] + t1.mat[1]*t2.mat[6] +
		       t1.mat[2]*t2.mat[10] + t1.mat[3]*t2.mat[14];
		tmp3 = t1.mat[0]*t2.mat[3] + t1.mat[1]*t2.mat[7] +
		       t1.mat[2]*t2.mat[11] + t1.mat[3]*t2.mat[15];
		tmp4 = t1.mat[4]*t2.mat[0] + t1.mat[5]*t2.mat[4] +
		       t1.mat[6]*t2.mat[8] + t1.mat[7]*t2.mat[12];
		tmp5 = t1.mat[4]*t2.mat[1] + t1.mat[5]*t2.mat[5] +
		       t1.mat[6]*t2.mat[9] + t1.mat[7]*t2.mat[13];
		tmp6 = t1.mat[4]*t2.mat[2] + t1.mat[5]*t2.mat[6] +
		       t1.mat[6]*t2.mat[10] + t1.mat[7]*t2.mat[14];
		tmp7 = t1.mat[4]*t2.mat[3] + t1.mat[5]*t2.mat[7] +
		       t1.mat[6]*t2.mat[11] + t1.mat[7]*t2.mat[15];
		tmp8 = t1.mat[8]*t2.mat[0] + t1.mat[9]*t2.mat[4] +
		       t1.mat[10]*t2.mat[8] + t1.mat[11]*t2.mat[12];
		tmp9 = t1.mat[8]*t2.mat[1] + t1.mat[9]*t2.mat[5] +
		       t1.mat[10]*t2.mat[9] + t1.mat[11]*t2.mat[13];
		tmp10 = t1.mat[8]*t2.mat[2] + t1.mat[9]*t2.mat[6] +
		        t1.mat[10]*t2.mat[10] + t1.mat[11]*t2.mat[14];
		tmp11 = t1.mat[8]*t2.mat[3] + t1.mat[9]*t2.mat[7] +
		        t1.mat[10]*t2.mat[11] + t1.mat[11]*t2.mat[15];

		if (t1.isAffine()) {
		    mat[12] = t2.mat[12];
		    mat[13] = t2.mat[13];
		    mat[14] = t2.mat[14];
		    mat[15] = t2.mat[15];
		} else {
		    double tmp12 = t1.mat[12]*t2.mat[0] + t1.mat[13]*t2.mat[4] +
		                   t1.mat[14]*t2.mat[8] + t1.mat[15]*t2.mat[12];
		    double tmp13 = t1.mat[12]*t2.mat[1] + t1.mat[13]*t2.mat[5] +
		                   t1.mat[14]*t2.mat[9] + t1.mat[15]*t2.mat[13];
		    double tmp14 = t1.mat[12]*t2.mat[2] + t1.mat[13]*t2.mat[6] +
			           t1.mat[14]*t2.mat[10] + t1.mat[15]*t2.mat[14];
		    double tmp15 = t1.mat[12]*t2.mat[3] + t1.mat[13]*t2.mat[7] +
		                   t1.mat[14]*t2.mat[11] + t1.mat[15]*t2.mat[15];
		    mat[12] = tmp12;
		    mat[13] = tmp13;
		    mat[14] = tmp14;
		    mat[15] = tmp15;
		}
	    }
	    mat[0] = tmp0;
	    mat[1] = tmp1;
	    mat[2] = tmp2;
	    mat[3] = tmp3;
	    mat[4] = tmp4;
	    mat[5] = tmp5;
	    mat[6] = tmp6;
	    mat[7] = tmp7;
	    mat[8] = tmp8;
	    mat[9] = tmp9;
	    mat[10] = tmp10;
	    mat[11] = tmp11;
	}


	if (((t1.dirtyBits & CONGRUENT_BIT) == 0) &&
	    ((t1.type & CONGRUENT) != 0) &&
	    ((t2.dirtyBits & CONGRUENT_BIT) == 0) &&
	    ((t2.type & CONGRUENT) != 0)) {
	    type = t1.type & t2.type;
	    dirtyBits = t1.dirtyBits | t2.dirtyBits | CLASSIFY_BIT |
		        ROTSCALESVD_DIRTY | RIGID_BIT;
	} else {
	    if (aff) {
		dirtyBits = ORTHO_BIT | CONGRUENT_BIT | RIGID_BIT |
                            CLASSIFY_BIT | ROTSCALESVD_DIRTY;
	    } else {
		dirtyBits = ALL_DIRTY;
	    }
	}

	if (autoNormalize) {
	    normalize();
	}
    }

    /**
     * Multiplies this transform by the inverse of transform t1. The final
     * value is placed into this matrix (this = this*t1^-1).
     * @param t1  the matrix whose inverse is computed.
     */
    public final void mulInverse(Transformation3D t1) {
	Transformation3D t2 = new Transformation3D();
	t2.autoNormalize = false;
	t2.invert(t1);
	this.mul(t2);
    }


    /**
     * Multiplies transform t1 by the inverse of transform t2. The final
     * value is placed into this matrix (this = t1*t2^-1).
     * @param t1  the left transform in the multiplication
     * @param t2  the transform whose inverse is computed.
     */
    public final void mulInverse(Transformation3D t1, Transformation3D t2) {
        Transformation3D t3 = new Transformation3D();
	t3.autoNormalize = false;
        t3.invert(t2);
        this.mul(t1,t3);
    }

    /**
     * Multiplies transform t1 by the transpose of transform t2 and places
     * the result into this transform (this = t1 * transpose(t2)).
     * @param t1  the transform on the left hand side of the multiplication
     * @param t2  the transform whose transpose is computed
     */
    public final void mulTransposeRight(Transformation3D t1, Transformation3D t2) {
	Transformation3D t3 = new Transformation3D();
	t3.autoNormalize = false;
	t3.transpose(t2);
	mul(t1, t3);
    }


    /**
     * Multiplies the transpose of transform t1 by transform t2 and places
     * the result into this matrix (this = transpose(t1) * t2).
     * @param t1  the transform whose transpose is computed
     * @param t2  the transform on the right hand side of the multiplication
     */
    public final void mulTransposeLeft(Transformation3D t1, Transformation3D t2){
	Transformation3D t3 = new Transformation3D();
	t3.autoNormalize = false;
	t3.transpose(t1);
	mul(t3, t2);
    }


    /**
     * Multiplies the transpose of transform t1 by the transpose of
     * transform t2 and places the result into this transform
     * (this = transpose(t1) * transpose(t2)).
     * @param t1  the transform on the left hand side of the multiplication
     * @param t2  the transform on the right hand side of the multiplication
     */
    public final void mulTransposeBoth(Transformation3D t1, Transformation3D t2) {
	Transformation3D t3 = new Transformation3D();
	Transformation3D t4 = new Transformation3D();
	t3.autoNormalize = false;
	t4.autoNormalize = false;
	t3.transpose(t1);
	t4.transpose(t2);
	mul(t3, t4);
    }


    /**
     * Normalizes the rotational components (upper 3x3) of this matrix
     * in place using a Singular Value Decomposition (SVD).
     * This operation ensures that the column vectors of this matrix
     * are orthogonal to each other.  The primary use of this method
     * is to correct for floating point errors that accumulate over
     * time when concatenating a large number of rotation matrices.
     * Note that the scale of the matrix is not altered by this method.
     */
    public final void normalize() {
        // Issue 253: Unable to normalize matrices with infinity or NaN
        if (!isAffine() && isInfOrNaN()) {
            return;
        }

	if ((dirtyBits & (ROTATION_BIT|SVD_BIT)) != 0) {
	    computeScaleRotation(true);
	} else 	if ((dirtyBits & (SCALE_BIT|SVD_BIT)) != 0) {
	    computeScales(true);
	}

	mat[0] = rot[0]*scales[0];
	mat[1] = rot[1]*scales[1];
	mat[2] = rot[2]*scales[2];
	mat[4] = rot[3]*scales[0];
	mat[5] = rot[4]*scales[1];
	mat[6] = rot[5]*scales[2];
	mat[8] = rot[6]*scales[0];
	mat[9] = rot[7]*scales[1];
	mat[10] = rot[8]*scales[2];
	dirtyBits |= CLASSIFY_BIT;
	dirtyBits &= ~ORTHO_BIT;
	type |= ORTHO;
    }

    /**
     * Normalizes the rotational components (upper 3x3) of transform t1
     * using a Singular Value Decomposition (SVD), and places the result
     * into this transform.
     * This operation ensures that the column vectors of this matrix
     * are orthogonal to each other.  The primary use of this method
     * is to correct for floating point errors that accumulate over
     * time when concatenating a large number of rotation matrices.
     * Note that the scale of the matrix is not altered by this method.
     *
     * @param t1  the source transform, which is not modified
     */
    public final void normalize(Transformation3D t1){
	set(t1);
	normalize();
    }

    /**
     * Sets the matrix, type, and state of this transform to the matrix,
     * type, and state of transform t1.
     * @param t1  the transform to be copied
     */
    public final void set(Transformation3D t1){
	mat[0] = t1.mat[0];
	mat[1] = t1.mat[1];
	mat[2] = t1.mat[2];
	mat[3] = t1.mat[3];
	mat[4] = t1.mat[4];
	mat[5] = t1.mat[5];
	mat[6] = t1.mat[6];
	mat[7] = t1.mat[7];
	mat[8] = t1.mat[8];
	mat[9] = t1.mat[9];
	mat[10] = t1.mat[10];
	mat[11] = t1.mat[11];
	mat[12] = t1.mat[12];
	mat[13] = t1.mat[13];
	mat[14] = t1.mat[14];
	mat[15] = t1.mat[15];
	type = t1.type;

	// don't copy rot[] and scales[]
	dirtyBits = t1.dirtyBits | ROTATION_BIT | SCALE_BIT;
        autoNormalize = t1.autoNormalize;
    }


    /**
     * Normalizes the rotational components (upper 3x3) of this transform
     * in place using a Cross Product (CP) normalization.
     * This operation ensures that the column vectors of this matrix
     * are orthogonal to each other.  The primary use of this method
     * is to correct for floating point errors that accumulate over
     * time when concatenating a large number of rotation matrices.
     * Note that the scale of the matrix is not altered by this method.
     */
    public final void normalizeCP()  {
        // Issue 253: Unable to normalize matrices with infinity or NaN
        if (!isAffine() && isInfOrNaN()) {
            return;
        }

	if ((dirtyBits & SCALE_BIT) != 0) {
	    computeScales(false);
	}

	double mag = mat[0]*mat[0] + mat[4]*mat[4] +
  	             mat[8]*mat[8];

	if (mag != 0) {
	    mag = 1.0/Math.sqrt(mag);
	    mat[0] = mat[0]*mag;
	    mat[4] = mat[4]*mag;
	    mat[8] = mat[8]*mag;
	}

	mag = mat[1]*mat[1] + mat[5]*mat[5] +
	      mat[9]*mat[9];

	if (mag != 0) {
	    mag = 1.0/Math.sqrt(mag);
	    mat[1] = mat[1]*mag;
	    mat[5] = mat[5]*mag;
	    mat[9] = mat[9]*mag;
	}
	mat[2] = (mat[4]*mat[9] - mat[5]*mat[8])*scales[0];
	mat[6] = (mat[1]*mat[8] - mat[0]*mat[9])*scales[1];
	mat[10] = (mat[0]*mat[5] - mat[1]*mat[4])*scales[2];

	mat[0] *= scales[0];
	mat[1] *= scales[0];
	mat[4] *= scales[1];
	mat[5] *= scales[1];
	mat[8] *= scales[2];
	mat[9] *= scales[2];

	// leave the AFFINE bit
	dirtyBits |= CONGRUENT_BIT | RIGID_BIT | CLASSIFY_BIT | ROTATION_BIT | SVD_BIT;
	dirtyBits &= ~ORTHO_BIT;
	type |= ORTHO;
    }


    /**
     * Normalizes the rotational components (upper 3x3) of transform t1
     * using a Cross Product (CP) normalization, and
     * places the result into this transform.
     * This operation ensures that the column vectors of this matrix
     * are orthogonal to each other.  The primary use of this method
     * is to correct for floating point errors that accumulate over
     * time when concatenating a large number of rotation matrices.
     * Note that the scale of the matrix is not altered by this method.
     *
     * @param t1 the transform to be normalized
     */
    public final void normalizeCP(Transformation3D t1) {
	set(t1);
	normalizeCP();
    }


    /**
     * Returns true if all of the data members of transform t1 are
     * equal to the corresponding data members in this Transformation3D.
     * @param t1  the transform with which the comparison is made
     * @return  true or false
     */
    public boolean equals(Transformation3D t1) {
	return (t1 != null) &&
	       (mat[0] == t1.mat[0]) && (mat[1] == t1.mat[1]) &&
	       (mat[2] == t1.mat[2]) && (mat[3] == t1.mat[3]) &&
	       (mat[4] == t1.mat[4]) && (mat[5] == t1.mat[5]) &&
	       (mat[6] == t1.mat[6]) && (mat[7] == t1.mat[7]) &&
	       (mat[8] == t1.mat[8]) && (mat[9] == t1.mat[9]) &&
	       (mat[10] == t1.mat[10]) && (mat[11] == t1.mat[11]) &&
	       (mat[12] == t1.mat[12]) && (mat[13] == t1.mat[13]) &&
	       (mat[14] == t1.mat[14]) && ( mat[15] == t1.mat[15]);
    }


   /**
     * Returns true if the Object o1 is of type Transformation3D and all of the
     * data members of o1 are equal to the corresponding data members in
     * this Transformation3D.
     * @param o1  the object with which the comparison is made.
     * @return  true or false
     */
    @Override
	public boolean equals(Object o1) {
	return (o1 instanceof Transformation3D) && equals((Transformation3D) o1);
    }


    /**
     * Returns true if the L-infinite distance between this matrix
     * and matrix m1 is less than or equal to the epsilon parameter,
     * otherwise returns false.  The L-infinite
     * distance is equal to
     * MAX[i=0,1,2,3 ; j=0,1,2,3 ; abs[(this.m(i,j) - m1.m(i,j)]
     * @param t1  the transform to be compared to this transform
     * @param epsilon  the threshold value
     */
    public boolean epsilonEquals(Transformation3D t1, double epsilon) {
        double diff;

        for (int i=0 ; i<16 ; i++) {
	    diff = mat[i] - t1.mat[i];
	    if ((diff < 0 ? -diff : diff) > epsilon) {
		return false;
	    }
        }
        return true;
    }


    /**
     * Transforms the normal parameter by this transform and places the value
     * back into normal.  The fourth element of the normal is assumed to be zero.
     * @param normal   the input normal to be transformed
     */
    public final void transform(IVecteur3d normal) {
        double x =  mat[0]*normal.obtenirX() + mat[1]*normal.obtenirY() + mat[2]*normal.obtenirZ();
        double y =  mat[4]*normal.obtenirX() + mat[5]*normal.obtenirY() + mat[6]*normal.obtenirZ();
        normal.définirZ(mat[8]*normal.obtenirX() + mat[9]*normal.obtenirY() + mat[10]*normal.obtenirZ());
        normal.définirX(x);
        normal.définirY(y);
    }

    /**
     *  Scales transform t1 by a Uniform scale matrix with scale
     *  factor s and then adds transform t2 (this = S*t1 + t2).
     *  @param s  the scale factor
     *  @param t1 the transform to be scaled
     *  @param t2 the transform to be added
     */
    public final void scaleAdd(double s, Transformation3D t1, Transformation3D t2) {
	for (int i=0 ; i<16 ; i++) {
	   mat[i] = s*t1.mat[i] + t2.mat[i];
	}

	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}
    }


    /**
     *  Scales this transform by a Uniform scale matrix with scale factor
     *  s and then adds transform t1 (this = S*this + t1).
     *  @param s  the scale factor
     *  @param t1 the transform to be added
     */
    public final void scaleAdd(double s, Transformation3D t1) {
	for (int i=0 ; i<16 ; i++) {
	    mat[i] = s*mat[i] + t1.mat[i];
	}

	dirtyBits = ALL_DIRTY;

	if (autoNormalize) {
	    normalize();
	}
   }


    /**
     * Creates an orthographic projection transform that mimics a standard,
     * camera-based,
     * view-model.  This transform maps coordinates from Eye Coordinates (EC)
     * to Clipping Coordinates (CC).  Note that unlike the similar function
     * in OpenGL, the clipping coordinates generated by the resulting
     * transform are in a right-handed coordinate system
     * (as are all other coordinate systems in Java 3D).
     * @param left the vertical line on the left edge of the near
     * clipping plane mapped to the left edge of the graphics window
     * @param right the vertical line on the right edge of the near
     * clipping plane mapped to the right edge of the graphics window
     * @param bottom the horizontal line on the bottom edge of the near
     * clipping plane mapped to the bottom edge of the graphics window
     * @param top the horizontal line on the top edge of the near
     * clipping plane mapped to the top edge of the graphics window
     * @param near the distance to the frustum's near clipping plane
     * (the value -near is the location of the near clip plane)
     * @param far the distance to the frustum's far clipping plane
     */
    public void ortho(double left, double right, double bottom,
                        double top, double near, double far) {
	double deltax = 1/(right - left);
	double deltay = 1/(top - bottom);
	double deltaz = 1/(far - near);

//	if ((deltax == 0.0) || (deltay == 0.0) || (deltaz == 0.0)) {
//	    return;
//	}

	mat[0] = 2.0 * deltax;
	mat[3] = -(right + left) * deltax;
	mat[5] = 2.0 * deltay;
	mat[7] = -(top + bottom) * deltay;
	mat[10] = 2.0 * deltaz;
	mat[11] = (far + near) * deltaz;
	mat[1] = mat[2] =  mat[4] = mat[6] = mat[8] =
	    mat[9] = mat[12] = mat[13] = mat[14] = 0;
	    mat[15] = 1;

	// Issue 253: set all dirty bits
	dirtyBits = ALL_DIRTY;
    }

    /**
     * get the scaling factor of matrix in this transform,
     * use for distance scaling
     */
    double getDistanceScale() {
	// The caller know that this matrix is affine
	// orthogonal before invoke this procedure

       if ((dirtyBits & SCALE_BIT) != 0) {
	   double max = mat[0]*mat[0] + mat[4]*mat[4] +
		        mat[8]*mat[8];
	   if (((dirtyBits & CONGRUENT_BIT) == 0) &&
	       ((type & CONGRUENT) != 0)) {
	       // in most case it is congruent
	       return Math.sqrt(max);
	   }
	   double tmp =  mat[1]*mat[1] + mat[5]*mat[5] +
		         mat[9]*mat[9];
	   if (tmp > max) {
	       max = tmp;
	   }
	   tmp = mat[2]*mat[2] + mat[6]*mat[6] + mat[10]*mat[10];
	   return Math.sqrt((tmp > max) ? tmp : max);
       }
       return max3(scales);
    }


    static private void  mat_mul(double[] m1, double[] m2, double[] m3) {

	double[] result = m3;
	if ((m1 == m3) || (m2 == m3)) {
	    result = new double[9];
	}

	result[0] =  m1[0]*m2[0] + m1[1]*m2[3] + m1[2]*m2[6];
	result[1] =  m1[0]*m2[1] + m1[1]*m2[4] + m1[2]*m2[7];
	result[2] =  m1[0]*m2[2] + m1[1]*m2[5] + m1[2]*m2[8];

	result[3] =  m1[3]*m2[0] + m1[4]*m2[3] + m1[5]*m2[6];
	result[4] =  m1[3]*m2[1] + m1[4]*m2[4] + m1[5]*m2[7];
	result[5] =  m1[3]*m2[2] + m1[4]*m2[5] + m1[5]*m2[8];

	result[6] =  m1[6]*m2[0] + m1[7]*m2[3] + m1[8]*m2[6];
	result[7] =  m1[6]*m2[1] + m1[7]*m2[4] + m1[8]*m2[7];
	result[8] =  m1[6]*m2[2] + m1[7]*m2[5] + m1[8]*m2[8];

	if (result != m3) {
	    for(int i=0;i<9;i++) {
		m3[i] = result[i];
	    }
	}
    }

    static private void  transpose_mat(double[] in, double[] out) {
	out[0] = in[0];
	out[1] = in[3];
	out[2] = in[6];

	out[3] = in[1];
	out[4] = in[4];
	out[5] = in[7];

	out[6] = in[2];
	out[7] = in[5];
	out[8] = in[8];
    }

    private void compute_svd(Transformation3D matrix, double[] outScale,
			     double[] outRot) {

	int i;
	double g;
	double m[] = new double[9];

	// if (!svdAllocd) {
	double[] u1 = new double[9];
	double[] v1 = new double[9];
	double[] t1 = new double[9];
	double[] t2 = new double[9];
	// double[] ts = new double[9];
	// double[] svdTmp = new double[9]; It is replaced by t1
	double[] svdRot = new double[9];
	// double[] single_values = new double[3]; replaced by t2

	double[] e = new double[3];
	double[] svdScales = new double[3];


	// XXXX: initialize to 0's if alread allocd? Should not have to, since
	// no operations depend on these being init'd to zero.

	int negCnt=0;
	double c1,c2,c3,c4;
	double s1,s2,s3,s4;


        svdRot[0] = m[0] = matrix.mat[0];
	svdRot[1] = m[1] = matrix.mat[1];
	svdRot[2] = m[2] = matrix.mat[2];
	svdRot[3] = m[3] = matrix.mat[4];
	svdRot[4] = m[4] = matrix.mat[5];
	svdRot[5] = m[5] = matrix.mat[6];
	svdRot[6] = m[6] = matrix.mat[8];
	svdRot[7] = m[7] = matrix.mat[9];
	svdRot[8] = m[8] = matrix.mat[10];

	// u1

	if( m[3]*m[3] < EPS ) {
	    u1[0] = 1.0; u1[1] = 0.0; u1[2] = 0.0;
	    u1[3] = 0.0; u1[4] = 1.0; u1[5] = 0.0;
	    u1[6] = 0.0; u1[7] = 0.0; u1[8] = 1.0;
	} else if( m[0]*m[0] < EPS ) {
	    t1[0] = m[0];
	    t1[1] = m[1];
	    t1[2] = m[2];
	    m[0] = m[3];
	    m[1] = m[4];
	    m[2] = m[5];

	    m[3] = -t1[0]; // zero
	    m[4] = -t1[1];
	    m[5] = -t1[2];

	    u1[0] =  0.0; u1[1] = 1.0;  u1[2] = 0.0;
	    u1[3] = -1.0; u1[4] = 0.0;  u1[5] = 0.0;
	    u1[6] =  0.0; u1[7] = 0.0;  u1[8] = 1.0;
	} else {
	    g = 1.0/Math.sqrt(m[0]*m[0] + m[3]*m[3]);
	    c1 = m[0]*g;
	    s1 = m[3]*g;
	    t1[0] = c1*m[0] + s1*m[3];
	    t1[1] = c1*m[1] + s1*m[4];
	    t1[2] = c1*m[2] + s1*m[5];

	    m[3] = -s1*m[0] + c1*m[3]; // zero
	    m[4] = -s1*m[1] + c1*m[4];
	    m[5] = -s1*m[2] + c1*m[5];

	    m[0] = t1[0];
	    m[1] = t1[1];
	    m[2] = t1[2];
	    u1[0] = c1;  u1[1] = s1;  u1[2] = 0.0;
	    u1[3] = -s1; u1[4] = c1;  u1[5] = 0.0;
	    u1[6] = 0.0; u1[7] = 0.0; u1[8] = 1.0;
	}

	// u2

	if( m[6]*m[6] < EPS  ) {
	} else if( m[0]*m[0] < EPS ){
	    t1[0] = m[0];
	    t1[1] = m[1];
	    t1[2] = m[2];
	    m[0] = m[6];
	    m[1] = m[7];
	    m[2] = m[8];

	    m[6] = -t1[0]; // zero
	    m[7] = -t1[1];
	    m[8] = -t1[2];

	    t1[0] = u1[0];
	    t1[1] = u1[1];
	    t1[2] = u1[2];
	    u1[0] = u1[6];
	    u1[1] = u1[7];
	    u1[2] = u1[8];

	    u1[6] = -t1[0]; // zero
	    u1[7] = -t1[1];
	    u1[8] = -t1[2];
	} else {
	    g = 1.0/Math.sqrt(m[0]*m[0] + m[6]*m[6]);
	    c2 = m[0]*g;
	    s2 = m[6]*g;
	    t1[0] = c2*m[0] + s2*m[6];
	    t1[1] = c2*m[1] + s2*m[7];
	    t1[2] = c2*m[2] + s2*m[8];

	    m[6] = -s2*m[0] + c2*m[6];
	    m[7] = -s2*m[1] + c2*m[7];
	    m[8] = -s2*m[2] + c2*m[8];
	    m[0] = t1[0];
	    m[1] = t1[1];
	    m[2] = t1[2];

	    t1[0] = c2*u1[0];
	    t1[1] = c2*u1[1];
	    u1[2]  = s2;

	    t1[6] = -u1[0]*s2;
	    t1[7] = -u1[1]*s2;
	    u1[8] = c2;
	    u1[0] = t1[0];
	    u1[1] = t1[1];
	    u1[6] = t1[6];
	    u1[7] = t1[7];
	}

	// v1

	if( m[2]*m[2] < EPS ) {
	    v1[0] = 1.0; v1[1] = 0.0; v1[2] = 0.0;
	    v1[3] = 0.0; v1[4] = 1.0; v1[5] = 0.0;
	    v1[6] = 0.0; v1[7] = 0.0; v1[8] = 1.0;
	} else if( m[1]*m[1] < EPS ) {
	    t1[2] = m[2];
	    t1[5] = m[5];
	    t1[8] = m[8];
	    m[2] = -m[1];
	    m[5] = -m[4];
	    m[8] = -m[7];

	    m[1] = t1[2]; // zero
	    m[4] = t1[5];
	    m[7] = t1[8];

	    v1[0] =  1.0; v1[1] = 0.0;  v1[2] = 0.0;
	    v1[3] =  0.0; v1[4] = 0.0;  v1[5] =-1.0;
	    v1[6] =  0.0; v1[7] = 1.0;  v1[8] = 0.0;
	} else {
	    g = 1.0/Math.sqrt(m[1]*m[1] + m[2]*m[2]);
	    c3 = m[1]*g;
	    s3 = m[2]*g;
	    t1[1] = c3*m[1] + s3*m[2];  // can assign to m[1]?
	    m[2] =-s3*m[1] + c3*m[2];  // zero
	    m[1] = t1[1];

	    t1[4] = c3*m[4] + s3*m[5];
	    m[5] =-s3*m[4] + c3*m[5];
	    m[4] = t1[4];

	    t1[7] = c3*m[7] + s3*m[8];
	    m[8] =-s3*m[7] + c3*m[8];
	    m[7] = t1[7];

	    v1[0] = 1.0; v1[1] = 0.0; v1[2] = 0.0;
	    v1[3] = 0.0; v1[4] =  c3; v1[5] = -s3;
	    v1[6] = 0.0; v1[7] =  s3; v1[8] =  c3;
	}

	// u3

	if( m[7]*m[7] < EPS ) {
	} else if( m[4]*m[4] < EPS ) {
	    t1[3] = m[3];
	    t1[4] = m[4];
	    t1[5] = m[5];
	    m[3] = m[6];   // zero
	    m[4] = m[7];
	    m[5] = m[8];

	    m[6] = -t1[3]; // zero
	    m[7] = -t1[4]; // zero
	    m[8] = -t1[5];

	    t1[3] = u1[3];
	    t1[4] = u1[4];
	    t1[5] = u1[5];
	    u1[3] = u1[6];
	    u1[4] = u1[7];
	    u1[5] = u1[8];

	    u1[6] = -t1[3]; // zero
	    u1[7] = -t1[4];
	    u1[8] = -t1[5];

	} else {
	    g = 1.0/Math.sqrt(m[4]*m[4] + m[7]*m[7]);
	    c4 = m[4]*g;
	    s4 = m[7]*g;
	    t1[3] = c4*m[3] + s4*m[6];
	    m[6] =-s4*m[3] + c4*m[6];  // zero
	    m[3] = t1[3];

	    t1[4] = c4*m[4] + s4*m[7];
	    m[7] =-s4*m[4] + c4*m[7];
	    m[4] = t1[4];

	    t1[5] = c4*m[5] + s4*m[8];
	    m[8] =-s4*m[5] + c4*m[8];
	    m[5] = t1[5];

	    t1[3] = c4*u1[3] + s4*u1[6];
	    u1[6] =-s4*u1[3] + c4*u1[6];
	    u1[3] = t1[3];

	    t1[4] = c4*u1[4] + s4*u1[7];
	    u1[7] =-s4*u1[4] + c4*u1[7];
	    u1[4] = t1[4];

	    t1[5] = c4*u1[5] + s4*u1[8];
	    u1[8] =-s4*u1[5] + c4*u1[8];
	    u1[5] = t1[5];
	}

	t2[0] = m[0];
	t2[1] = m[4];
	t2[2] = m[8];
	e[0] = m[1];
	e[1] = m[5];

	if( e[0]*e[0]>EPS || e[1]*e[1]>EPS ) {
	    compute_qr( t2, e, u1, v1);
	}

	svdScales[0] = t2[0];
	svdScales[1] = t2[1];
	svdScales[2] = t2[2];


	// Do some optimization here. If scale is unity, simply return the rotation matric.
	if(almostOne(Math.abs(svdScales[0])) &&
	   almostOne(Math.abs(svdScales[1])) &&
	   almostOne(Math.abs(svdScales[2]))) {

	    for(i=0;i<3;i++)
		if(svdScales[i]<0.0)
		    negCnt++;

	    if((negCnt==0)||(negCnt==2)) {
		//System.err.println("Optimize!!");
		outScale[0] = outScale[1] = outScale[2] = 1.0;
		for(i=0;i<9;i++)
		    outRot[i] = svdRot[i];

		return;
	    }
	}

	// XXXX: could eliminate use of t1 and t1 by making a new method which
	// transposes and multiplies two matricies
	transpose_mat(u1, t1);
	transpose_mat(v1, t2);


	svdReorder( m, t1, t2, svdRot, svdScales, outRot, outScale);
    }


    private void svdReorder( double[] m, double[] t1, double[] t2, double[] rot,
			     double[] scales, double[] outRot, double[] outScale) {

	int in0, in1, in2, index,i;
	int[] svdOut = new int[3];
	double[] svdMag = new double[3];


	// check for rotation information in the scales
	if(scales[0] < 0.0 ) {   // move the rotation info to rotation matrix
	    scales[0] = -scales[0];
	    t2[0] = -t2[0];
	    t2[1] = -t2[1];
	    t2[2] = -t2[2];
	}
	if(scales[1] < 0.0 ) {   // move the rotation info to rotation matrix
	    scales[1] = -scales[1];
	    t2[3] = -t2[3];
	    t2[4] = -t2[4];
	    t2[5] = -t2[5];
	}
	if(scales[2] < 0.0 ) {   // move the rotation info to rotation matrix
	    scales[2] = -scales[2];
	    t2[6] = -t2[6];
	    t2[7] = -t2[7];
	    t2[8] = -t2[8];
	}


	mat_mul(t1,t2,rot);

	// check for equal scales case  and do not reorder
	if(almostEqual(Math.abs(scales[0]), Math.abs(scales[1])) &&
	   almostEqual(Math.abs(scales[1]), Math.abs(scales[2]))   ){
	    for(i=0;i<9;i++){
		outRot[i] = rot[i];
	    }
	    for(i=0;i<3;i++){
		outScale[i] = scales[i];
	    }

	}else {

	    // sort the order of the results of SVD
	    if( scales[0] > scales[1]) {
		if( scales[0] > scales[2] ) {
		    if( scales[2] > scales[1] ) {
			svdOut[0] = 0; svdOut[1] = 2; svdOut[2] = 1; // xzy
		    } else {
			svdOut[0] = 0; svdOut[1] = 1; svdOut[2] = 2; // xyz
		    }
		} else {
		    svdOut[0] = 2; svdOut[1] = 0; svdOut[2] = 1; // zxy
		}
	    } else {  // y > x
		if( scales[1] > scales[2] ) {
		    if( scales[2] > scales[0] ) {
			svdOut[0] = 1; svdOut[1] = 2; svdOut[2] = 0; // yzx
		    } else {
			svdOut[0] = 1; svdOut[1] = 0; svdOut[2] = 2; // yxz
		    }
		} else  {
		    svdOut[0] = 2; svdOut[1] = 1; svdOut[2] = 0; // zyx
		}
	    }


	    // sort the order of the input matrix
	    svdMag[0] = (m[0]*m[0] + m[1]*m[1] + m[2]*m[2]);
	    svdMag[1] = (m[3]*m[3] + m[4]*m[4] + m[5]*m[5]);
	    svdMag[2] = (m[6]*m[6] + m[7]*m[7] + m[8]*m[8]);


	    if( svdMag[0] > svdMag[1]) {
		if( svdMag[0] > svdMag[2] ) {
		    if( svdMag[2] > svdMag[1] )  {
			// 0 - 2 - 1
			in0 = 0; in2 = 1; in1 = 2;// xzy
		    } else {
			// 0 - 1 - 2
			in0 = 0; in1 = 1; in2 = 2; // xyz
		    }
		} else {
		    // 2 - 0 - 1
		    in2 = 0; in0 = 1; in1 = 2;  // zxy
		}
	    } else {  // y > x   1>0
		if( svdMag[1] > svdMag[2] ) {  // 1>2
		    if( svdMag[2] > svdMag[0] )  { // 2>0
			// 1 - 2 - 0
			in1 = 0; in2 = 1; in0 = 2; // yzx
		    } else {
			// 1 - 0 - 2
			in1 = 0; in0 = 1; in2 = 2; // yxz
		    }
		} else  {
		    // 2 - 1 - 0
		    in2 = 0; in1 = 1; in0 = 2; // zyx
		}
	    }


	    index = svdOut[in0];
	    outScale[0] = scales[index];

	    index = svdOut[in1];
	    outScale[1] = scales[index];

	    index = svdOut[in2];
	    outScale[2] = scales[index];

	    index = svdOut[in0];

	    outRot[0] = rot[index];

	    index = svdOut[in0]+3;
	    outRot[0+3] = rot[index];

	    index = svdOut[in0]+6;
	    outRot[0+6] = rot[index];

	    index = svdOut[in1];
	    outRot[1] = rot[index];

	    index = svdOut[in1]+3;
	    outRot[1+3] = rot[index];

	    index = svdOut[in1]+6;
	    outRot[1+6] = rot[index];

	    index = svdOut[in2];
	    outRot[2] = rot[index];

	    index = svdOut[in2]+3;
	    outRot[2+3] = rot[index];

	    index = svdOut[in2]+6;
	    outRot[2+6] = rot[index];
	}

    }

    private int compute_qr( double[] s, double[] e, double[] u, double[] v) {
	int k;
	boolean converged;
	double shift,r;

	double utemp,vtemp;
	double f,g;

	final int MAX_INTERATIONS = 10;
	final double CONVERGE_TOL = 4.89E-15;

	double[]   cosl  = new double[2];
	double[]   cosr  = new double[2];
	double[]   sinl  = new double[2];
	double[]   sinr  = new double[2];
	double[]   qr_m  = new double[9];


	double c_b48 = 1.;
	int first;
	converged = false;

	first = 1;

	if( Math.abs(e[1]) < CONVERGE_TOL || Math.abs(e[0]) < CONVERGE_TOL) converged = true;

	for(k=0;k<MAX_INTERATIONS && !converged;k++) {
	    shift = compute_shift( s[1], e[1], s[2]);
	    f = (Math.abs(s[0]) - shift) * (d_sign(c_b48, s[0]) + shift/s[0]);
	    g = e[0];
	    r = compute_rot(f, g, sinr, cosr,  0, first);
	    f = cosr[0] * s[0] + sinr[0] * e[0];
	    e[0] = cosr[0] * e[0] - sinr[0] * s[0];
	    g = sinr[0] * s[1];
	    s[1] = cosr[0] * s[1];

	    r = compute_rot(f, g, sinl, cosl, 0, first);
	    first = 0;
	    s[0] = r;
	    f = cosl[0] * e[0] + sinl[0] * s[1];
	    s[1] = cosl[0] * s[1] - sinl[0] * e[0];
	    g = sinl[0] * e[1];
	    e[1] =  cosl[0] * e[1];

	    r = compute_rot(f, g, sinr, cosr, 1, first);
	    e[0] = r;
	    f = cosr[1] * s[1] + sinr[1] * e[1];
	    e[1] = cosr[1] * e[1] - sinr[1] * s[1];
	    g = sinr[1] * s[2];
	    s[2] = cosr[1] * s[2];

	    r = compute_rot(f, g, sinl, cosl, 1, first);
	    s[1] = r;
	    f = cosl[1] * e[1] + sinl[1] * s[2];
	    s[2] = cosl[1] * s[2] - sinl[1] * e[1];
	    e[1] = f;

	    // update u  matrices
	    utemp = u[0];
	    u[0] = cosl[0]*utemp + sinl[0]*u[3];
	    u[3] = -sinl[0]*utemp + cosl[0]*u[3];
	    utemp = u[1];
	    u[1] = cosl[0]*utemp + sinl[0]*u[4];
	    u[4] = -sinl[0]*utemp + cosl[0]*u[4];
	    utemp = u[2];
	    u[2] = cosl[0]*utemp + sinl[0]*u[5];
	    u[5] = -sinl[0]*utemp + cosl[0]*u[5];

	    utemp = u[3];
	    u[3] = cosl[1]*utemp + sinl[1]*u[6];
	    u[6] = -sinl[1]*utemp + cosl[1]*u[6];
	    utemp = u[4];
	    u[4] = cosl[1]*utemp + sinl[1]*u[7];
	    u[7] = -sinl[1]*utemp + cosl[1]*u[7];
	    utemp = u[5];
	    u[5] = cosl[1]*utemp + sinl[1]*u[8];
	    u[8] = -sinl[1]*utemp + cosl[1]*u[8];

	    // update v  matrices

	    vtemp = v[0];
	    v[0] = cosr[0]*vtemp + sinr[0]*v[1];
	    v[1] = -sinr[0]*vtemp + cosr[0]*v[1];
	    vtemp = v[3];
	    v[3] = cosr[0]*vtemp + sinr[0]*v[4];
	    v[4] = -sinr[0]*vtemp + cosr[0]*v[4];
	    vtemp = v[6];
	    v[6] = cosr[0]*vtemp + sinr[0]*v[7];
	    v[7] = -sinr[0]*vtemp + cosr[0]*v[7];

	    vtemp = v[1];
	    v[1] = cosr[1]*vtemp + sinr[1]*v[2];
	    v[2] = -sinr[1]*vtemp + cosr[1]*v[2];
	    vtemp = v[4];
	    v[4] = cosr[1]*vtemp + sinr[1]*v[5];
	    v[5] = -sinr[1]*vtemp + cosr[1]*v[5];
	    vtemp = v[7];
	    v[7] = cosr[1]*vtemp + sinr[1]*v[8];
	    v[8] = -sinr[1]*vtemp + cosr[1]*v[8];

	    // if(debug)System.err.println("\n*********************** iteration #"+k+" ***********************\n");

	    qr_m[0] = s[0];  qr_m[1] = e[0]; qr_m[2] = 0.0;
	    qr_m[3] =  0.0;  qr_m[4] = s[1]; qr_m[5] =e[1];
	    qr_m[6] =  0.0;  qr_m[7] =  0.0; qr_m[8] =s[2];

	    if( Math.abs(e[1]) < CONVERGE_TOL || Math.abs(e[0]) < CONVERGE_TOL) converged = true;
	}

	if( Math.abs(e[1]) < CONVERGE_TOL ) {
	    compute_2X2( s[0],e[0],s[1],s,sinl,cosl,sinr,cosr, 0);

	    utemp = u[0];
	    u[0] = cosl[0]*utemp + sinl[0]*u[3];
	    u[3] = -sinl[0]*utemp + cosl[0]*u[3];
	    utemp = u[1];
	    u[1] = cosl[0]*utemp + sinl[0]*u[4];
	    u[4] = -sinl[0]*utemp + cosl[0]*u[4];
	    utemp = u[2];
	    u[2] = cosl[0]*utemp + sinl[0]*u[5];
	    u[5] = -sinl[0]*utemp + cosl[0]*u[5];

	    // update v  matrices

	    vtemp = v[0];
	    v[0] = cosr[0]*vtemp + sinr[0]*v[1];
	    v[1] = -sinr[0]*vtemp + cosr[0]*v[1];
	    vtemp = v[3];
	    v[3] = cosr[0]*vtemp + sinr[0]*v[4];
	    v[4] = -sinr[0]*vtemp + cosr[0]*v[4];
	    vtemp = v[6];
	    v[6] = cosr[0]*vtemp + sinr[0]*v[7];
	    v[7] = -sinr[0]*vtemp + cosr[0]*v[7];
	} else {
	    compute_2X2( s[1],e[1],s[2],s,sinl,cosl,sinr,cosr,1);

	    utemp = u[3];
	    u[3] = cosl[0]*utemp + sinl[0]*u[6];
	    u[6] = -sinl[0]*utemp + cosl[0]*u[6];
	    utemp = u[4];
	    u[4] = cosl[0]*utemp + sinl[0]*u[7];
	    u[7] = -sinl[0]*utemp + cosl[0]*u[7];
	    utemp = u[5];
	    u[5] = cosl[0]*utemp + sinl[0]*u[8];
	    u[8] = -sinl[0]*utemp + cosl[0]*u[8];

	    // update v  matrices

	    vtemp = v[1];
	    v[1] = cosr[0]*vtemp + sinr[0]*v[2];
	    v[2] = -sinr[0]*vtemp + cosr[0]*v[2];
	    vtemp = v[4];
	    v[4] = cosr[0]*vtemp + sinr[0]*v[5];
	    v[5] = -sinr[0]*vtemp + cosr[0]*v[5];
	    vtemp = v[7];
	    v[7] = cosr[0]*vtemp + sinr[0]*v[8];
	    v[8] = -sinr[0]*vtemp + cosr[0]*v[8];
	}

	return(0);
    }

    static final double max( double a, double b) {
	return ( a > b ? a : b);
    }

    static final double min( double a, double b) {
	return ( a < b ? a : b);
    }

    static final double d_sign(double a, double b) {
        double x =  (a >= 0 ? a : - a);
	return( b >= 0 ? x : -x);
    }

    static final double compute_shift( double f, double g, double h) {
	double d__1, d__2;
	double fhmn, fhmx, c, fa, ga, ha, as, at, au;
	double ssmin;

	fa = Math.abs(f);
	ga = Math.abs(g);
	ha = Math.abs(h);
	fhmn = min(fa,ha);
	fhmx = max(fa,ha);
	if (fhmn == 0.) {
	    ssmin = 0.;
	    if (fhmx == 0.) {
	    } else {
		d__1 = min(fhmx,ga) / max(fhmx,ga);
	    }
	} else {
	    if (ga < fhmx) {
		as = fhmn / fhmx + 1.;
		at = (fhmx - fhmn) / fhmx;
		d__1 = ga / fhmx;
		au = d__1 * d__1;
		c = 2. / (Math.sqrt(as * as + au) + Math.sqrt(at * at + au));
		ssmin = fhmn * c;
	    } else {
		au = fhmx / ga;
		if (au == 0.) {


		    ssmin = fhmn * fhmx / ga;
		} else {
		    as = fhmn / fhmx + 1.;
		    at = (fhmx - fhmn) / fhmx;
		    d__1 = as * au;
		    d__2 = at * au;
		    c = 1. / (Math.sqrt(d__1 * d__1 + 1.) + Math.sqrt(d__2 * d__2 + 1.));
		    ssmin = fhmn * c * au;
		    ssmin += ssmin;
		}
	    }
	}

	return(ssmin);
    }

    static int compute_2X2( double f, double g, double h, double[] single_values,
			    double[] snl, double[] csl, double[] snr, double[] csr, int index)  {

	double c_b3 = 2.;
	double c_b4 = 1.;

	double d__1;
	int pmax;
	double temp;
	boolean swap;
	double a, d, l, m, r, s, t, tsign, fa, ga, ha;
	double ft, gt, ht, mm;
	boolean gasmal;
	double tt, clt, crt, slt, srt;
	double ssmin,ssmax;

	ssmax = single_values[0];
	ssmin = single_values[1];
	clt = 0.0;
	crt = 0.0;
	slt = 0.0;
	srt = 0.0;
	tsign = 0.0;

	ft = f;
	fa = Math.abs(ft);
	ht = h;
	ha = Math.abs(h);

	pmax = 1;
	if( ha > fa)
	    swap = true;
	else
	    swap = false;

	if (swap) {
	    pmax = 3;
	    temp = ft;
	    ft = ht;
	    ht = temp;
	    temp = fa;
	    fa = ha;
	    ha = temp;

	}
	gt = g;
	ga = Math.abs(gt);
	if (ga == 0.) {

	    single_values[1] = ha;
	    single_values[0] = fa;
	    clt = 1.;
	    crt = 1.;
	    slt = 0.;
	    srt = 0.;
	} else {
	    gasmal = true;

	    if (ga > fa) {
		pmax = 2;
		if (fa / ga < EPS) {

		    gasmal = false;
		    ssmax = ga;
		    if (ha > 1.) {
			ssmin = fa / (ga / ha);
		    } else {
			ssmin = fa / ga * ha;
		    }
		    clt = 1.;
		    slt = ht / gt;
		    srt = 1.;
		    crt = ft / gt;
		}
	    }
	    if (gasmal) {

		d = fa - ha;
		if (d == fa) {

		    l = 1.;
		} else {
		    l = d / fa;
		}

		m = gt / ft;

		t = 2. - l;

		mm = m * m;
		tt = t * t;
		s = Math.sqrt(tt + mm);

		if (l == 0.) {
		    r = Math.abs(m);
		} else {
		    r = Math.sqrt(l * l + mm);
		}

		a = (s + r) * .5;

		if (ga > fa) {
		    pmax = 2;
		    if (fa / ga < EPS) {

			gasmal = false;
			ssmax = ga;
			if (ha > 1.) {
			    ssmin = fa / (ga / ha);
			} else {
			    ssmin = fa / ga * ha;
			}
			clt = 1.;
			slt = ht / gt;
			srt = 1.;
			crt = ft / gt;
		    }
		}
		if (gasmal) {

		    d = fa - ha;
		    if (d == fa) {

			l = 1.;
		    } else {
			l = d / fa;
		    }

		    m = gt / ft;

		    t = 2. - l;

		    mm = m * m;
		    tt = t * t;
		    s = Math.sqrt(tt + mm);

		    if (l == 0.) {
			r = Math.abs(m);
		    } else {
			r = Math.sqrt(l * l + mm);
		    }

		    a = (s + r) * .5;


		    ssmin = ha / a;
		    ssmax = fa * a;
		    if (mm == 0.) {

			if (l == 0.) {
			    t = d_sign(c_b3, ft) * d_sign(c_b4, gt);
			} else {
			    t = gt / d_sign(d, ft) + m / t;
			}
		    } else {
			t = (m / (s + t) + m / (r + l)) * (a + 1.);
		    }
		    l = Math.sqrt(t * t + 4.);
		    crt = 2. / l;
		    srt = t / l;
		    clt = (crt + srt * m) / a;
		    slt = ht / ft * srt / a;
		}
	    }
	    if (swap) {
		csl[0] = srt;
		snl[0] = crt;
		csr[0] = slt;
		snr[0] = clt;
	    } else {
		csl[0] = clt;
		snl[0] = slt;
		csr[0] = crt;
		snr[0] = srt;
	    }

	    if (pmax == 1) {
		tsign = d_sign(c_b4, csr[0]) * d_sign(c_b4, csl[0]) * d_sign(c_b4, f);
	    }
	    if (pmax == 2) {
		tsign = d_sign(c_b4, snr[0]) * d_sign(c_b4, csl[0]) * d_sign(c_b4, g);
	    }
	    if (pmax == 3) {
		tsign = d_sign(c_b4, snr[0]) * d_sign(c_b4, snl[0]) * d_sign(c_b4, h);
	    }
	    single_values[index] = d_sign(ssmax, tsign);
	    d__1 = tsign * d_sign(c_b4, f) * d_sign(c_b4, h);
	    single_values[index+1] = d_sign(ssmin, d__1);


	}
	return 0;
    }

    static  double compute_rot( double f, double g, double[] sin, double[] cos, int index, int first) {
	double cs,sn;
	int i;
	double scale;
	int count;
	double f1, g1;
	double r;
	final double safmn2 = 2.002083095183101E-146;
	final double safmx2 = 4.994797680505588E+145;

	if (g == 0.) {
	    cs = 1.;
	    sn = 0.;
	    r = f;
	} else if (f == 0.) {
	    cs = 0.;
	    sn = 1.;
	    r = g;
	} else {
	    f1 = f;
	    g1 = g;
	    scale = max(Math.abs(f1),Math.abs(g1));
	    if (scale >= safmx2) {
		count = 0;
		while(scale >= safmx2) {
		    ++count;
		    f1 *= safmn2;
		    g1 *= safmn2;
		    scale = max(Math.abs(f1),Math.abs(g1));
		}
		r = Math.sqrt(f1*f1 + g1*g1);
		cs = f1 / r;
		sn = g1 / r;
		for (i = 1; i <= count; ++i) {
		    r *= safmx2;
		}
	    } else if (scale <= safmn2) {
		count = 0;
		while(scale <= safmn2) {
		    ++count;
		    f1 *= safmx2;
		    g1 *= safmx2;
		    scale = max(Math.abs(f1),Math.abs(g1));
		}
		r = Math.sqrt(f1*f1 + g1*g1);
		cs = f1 / r;
		sn = g1 / r;
		for (i = 1; i <= count; ++i) {
		    r *= safmn2;
		}
	    } else {
		r = Math.sqrt(f1*f1 + g1*g1);
		cs = f1 / r;
		sn = g1 / r;
	    }
	    if (Math.abs(f) > Math.abs(g) && cs < 0.) {
		cs = -cs;
		sn = -sn;
		r = -r;
	    }
	}
	sin[index] = sn;
	cos[index] = cs;
	return r;

    }

    static final private double max3( double[] values) {
	if( values[0] > values[1] ) {
	    if( values[0] > values[2] )
		return(values[0]);
	    else
		return(values[2]);
	} else {
	    if( values[1] > values[2] )
		return(values[1]);
	    else
		return(values[2]);
	}
    }


    final private void computeScales(boolean forceSVD) {

	if(scales == null)
	    scales = new double[3];

	if ((!forceSVD || ((dirtyBits & SVD_BIT) == 0)) && isAffine()) {
	    if (isCongruent()) {
		if (((dirtyBits & RIGID_BIT) == 0) &&
                    ((type & RIGID) != 0)) {
		    scales[0] = scales[1] = scales[2] = 1;
		    dirtyBits &= ~SCALE_BIT;
		    return;
		}
		scales[0] = scales[1] = scales[2] =
		    Math.sqrt(mat[0]*mat[0] + mat[4]*mat[4] +
			      mat[8]*mat[8]);
		dirtyBits &= ~SCALE_BIT;
		return;
	    }
	    if (isOrtho()) {
		scales[0] = Math.sqrt(mat[0]*mat[0] + mat[4]*mat[4] +
				      mat[8]*mat[8]);
		scales[1] = Math.sqrt(mat[1]*mat[1] + mat[5]*mat[5] +
				      mat[9]*mat[9]);
		scales[2] = Math.sqrt(mat[2]*mat[2] + mat[6]*mat[6] +
				      mat[10]*mat[10]);
		dirtyBits &= ~SCALE_BIT;
		return;
	    }
	}
	// fall back to use SVD decomposition
	if (rot == null)
	    rot = new double[9];

	compute_svd(this, scales, rot);
	dirtyBits &= ~ROTSCALESVD_DIRTY;
    }

    final private void computeScaleRotation(boolean forceSVD) {

	if(rot == null)
	    rot = new double[9];

	if(scales == null)
	    scales = new double[3];

	if ((!forceSVD || ((dirtyBits & SVD_BIT) == 0)) && isAffine()) {
	    if (isCongruent()) {
		if (((dirtyBits & RIGID_BIT) == 0) &&
                    ((type & RIGID) != 0)) {
		    rot[0] = mat[0];
		    rot[1] = mat[1];
		    rot[2] = mat[2];
		    rot[3] = mat[4];
		    rot[4] = mat[5];
		    rot[5] = mat[6];
		    rot[6] = mat[8];
		    rot[7] = mat[9];
		    rot[8] = mat[10];
		    scales[0] = scales[1] = scales[2] = 1;
		    dirtyBits &= (~ROTATION_BIT | ~SCALE_BIT);
		    return;
		}
		double s = Math.sqrt(mat[0]*mat[0] + mat[4]*mat[4] + mat[8]*mat[8]);
		if (s == 0) {
		    compute_svd(this, scales, rot);
		    return;
		}
		scales[0] = scales[1] = scales[2] = s;
		s = 1/s;
		rot[0] = mat[0]*s;
		rot[1] = mat[1]*s;
		rot[2] = mat[2]*s;
		rot[3] = mat[4]*s;
		rot[4] = mat[5]*s;
		rot[5] = mat[6]*s;
		rot[6] = mat[8]*s;
		rot[7] = mat[9]*s;
		rot[8] = mat[10]*s;
		dirtyBits &= (~ROTATION_BIT | ~SCALE_BIT);
		return;
	    }
	    if (isOrtho()) {
		double s;

		scales[0] = Math.sqrt(mat[0]*mat[0] + mat[4]*mat[4] + mat[8]*mat[8]);
		scales[1] = Math.sqrt(mat[1]*mat[1] + mat[5]*mat[5] + mat[9]*mat[9]);
		scales[2] = Math.sqrt(mat[2]*mat[2] + mat[6]*mat[6] + mat[10]*mat[10]);

		if ((scales[0] == 0) || (scales[1] == 0) || (scales[2] == 0)) {
		    compute_svd(this, scales, rot);
		    return;
		}
		s = 1/scales[0];
		rot[0] = mat[0]*s;
		rot[3] = mat[4]*s;
		rot[6] = mat[8]*s;
		s = 1/scales[1];
		rot[1] = mat[1]*s;
		rot[4] = mat[5]*s;
		rot[7] = mat[9]*s;
		s = 1/scales[2];
		rot[2] = mat[2]*s;
		rot[5] = mat[6]*s;
		rot[8] = mat[10]*s;
		dirtyBits &= (~ROTATION_BIT | ~SCALE_BIT);
		return;
	    }
	}
	// fall back to use SVD decomposition
	compute_svd(this, scales, rot);
	dirtyBits &= ~ROTSCALESVD_DIRTY;
    }


    final void getRotation(Transformation3D t) {
	if ((dirtyBits & ROTATION_BIT)!= 0) {
	    computeScaleRotation(false);
	}

        t.mat[3] = t.mat[7] = t.mat[11] = t.mat[12] = t.mat[13] =
	    t.mat[14] = 0;
	t.mat[15] = 1;
	t.mat[0] = rot[0];
	t.mat[1] = rot[1];
	t.mat[2] = rot[2];
	t.mat[4] = rot[3];
	t.mat[5] = rot[4];
	t.mat[6] = rot[5];
	t.mat[8] = rot[6];
	t.mat[9] = rot[7];
	t.mat[10] = rot[8];

	// Issue 253: set all dirty bits
	t.dirtyBits = ALL_DIRTY;
    }

    // somehow CanvasViewCache will directly modify mat[]
    // instead of calling ortho(). So we need to reset dirty bit
    final void setOrthoDirtyBit() {
	// Issue 253: set all dirty bits
	dirtyBits = ALL_DIRTY;
	type = 0;
    }

    // Fix for Issue 167 -- don't classify matrices with Infinity or NaN values
    // as affine
    private final boolean isInfOrNaN() {
        // The following is a faster version of the check.
        // Instead of 3 tests per array element (Double.isInfinite is 2 tests),
        // for a total of 48 tests, we will do 16 multiplies and 1 test.
        double d = 0.0;
        for (int i = 0; i < 16; i++) {
            d *= mat[i];
        }

        return d != 0.0;
    }

    private boolean isInfOrNaN(double val) {
        return Double.isNaN(val) || Double.isInfinite(val);
    }

    private boolean isInfOrNaN(IVecteur3d v) {
        return (Double.isNaN(v.obtenirX()) || Double.isInfinite(v.obtenirX()) ||
        		Double.isNaN(v.obtenirY()) || Double.isInfinite(v.obtenirY()) ||
        		Double.isNaN(v.obtenirZ()) || Double.isInfinite(v.obtenirZ()));
    }
    
}// FIN Classe
